\documentclass{article}
\usepackage{fullpage}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{courier}

% If your system does not have the AMS fonts version 2.0 installed, then
% remove the useAMS option.
%
% useAMS allows you to obtain upright Greek characters.
% e.g. \umu, \upi etc.  See the section on "Upright Greek characters" in
% this guide for further information.
%
% If you are using AMS 2.0 fonts, bold math letters/symbols are available
% at a larger range of sizes for NFSS release 1 and 2 (using \boldmath or
% preferably \bmath).
%
% The usenatbib command allows the use of Patrick Daly's natbib.sty for
% cross-referencing.
%
% If you wish to typeset the paper in Times font (if you do not have the
% PostScript Type 1 Computer Modern fonts you will need to do this to get
% smoother fonts in a PDF file) then uncomment the next line
% \usepackage{Times}

%%%%% AUTHORS - PLACE YOUR OWN MACROS HERE %%%%%


%% cosmology
%% ions / absorbers
\newcommand{\HI}{\hbox{{\sc H}{\sc i}} }
\newcommand{\NHI}{{N_{\rm HI}}}
\newcommand{\fNHI}{f(N_{\rm HI},X)}

\newcommand{\DXobs}{\Delta X^{\rm obs}}
\newcommand{\DXsim}{\Delta X^{\rm sim}}

\newcommand{\dXobs}{dX^{\rm obs}}
\newcommand{\dXsim}{dX^{\rm sim}}

\newcommand{\nel}{n_{\rm e} }
\newcommand{\nH}{n_{\rm _{H}} }
\newcommand{\nHzero}{n_{\rm _{H,0}} }
\newcommand{\nHI}{n_{\rm _{HI}} }
\newcommand{\nHII}{n_{\rm _{HII}} }
\newcommand{\xHI}{x_{\rm _{HI}} }
\newcommand{\xHIzero}{x_{\rm _{HI,0}} }
\newcommand{\xHII}{x_{\rm _{HII}} }
\newcommand{\xHIIzero}{x_{\rm _{HII,0}} }
\newcommand{\xo}{x_{\rm _{0}} }

\newcommand{\sigth}{\sigma_{\rm 0}}
\newcommand{\Inu}{I_{\nu} }      % Specific Intensity
\newcommand{\nuth}{\nu_{\rm th}}
\newcommand{\numax}{\nu_{\rm max}}
\newcommand{\Jnu}{J_{\nu}}
\newcommand{\signu}{\sigma_{\nu}}
\newcommand{\taunu}{\tau_{\nu}}

\newcommand{\Ob}{\Omega_{\rm b} }
\newcommand{\Om}{\Omega_{\rm m} }
\newcommand{\OHI}{\Omega_{\rm HI} }
\newcommand{\OHH}{\Omega_{\rm H2} }
\newcommand{\Ol}{\Omega_{\Lambda} }
\newcommand{\Mpch}{h^{-1} \rm{\,Mpc} }

% codes / simulations
\newcommand{\owls}{{\small OWLS} }
\newcommand{\urchin}{{\small URCHIN} }

% journals

\newcommand{\mnras}{MNRAS}
\newcommand{\apjl}{ApJL}
\newcommand{\apj}{ApJ}
\newcommand{\apjs}{ApJS}
\newcommand{\aap}{A\&A}
\newcommand{\araa}{ARA\&A}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Hydrogen Ionization States \\ Analytic Solutions}
\author{Gabriel Altay}
\begin{document}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

These notes are meant to document a series of analytic solutions for the evolution of the ionization state of hydrogen.  They consider a simple model atom that can only be in the ground state or ionized. 

\subsection{Basic Math}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In computational physics, we often deal with tabulated functions.  Commonly, the independent variable is sampled either linearly or logarithmically. Here we derive some simple relationships useful when doing integrals.  To be more concrete, suppose $w = \log x$ where $\log$ is the base 10 logarithm.  We then have $w = \log x = \ln x / \ln 10$ or $10^w = x$.  Now, suppose we want to know the following integral, $\int_a^b f(x) dx$ but the value of $f(x)$ is tabulated with logarithmic samples in $x$.  The function $f$ can always be evaluated either in terms of $x$ or $w$. However, to perform the integral we need to know the relationship between $dx$ and $dw$.  We can calculate this using the following facts,  
\begin{equation}
w \equiv \log x = \frac{\ln x}{\ln 10}, \quad
10^w = x, \quad
\frac{d}{dx} \ln x = \frac{1}{x}
\end{equation} 
Then the differentials are related as $dx = 10^w \ln 10 \, dw$ and the integral can be rewriten as, 
\begin{equation}
\int_a^b f(x) dx = 
\int_{\log a}^{\log b} f(w) 10^w \ln 10 \, dw
\end{equation} 



\subsection{Special Functions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Many of the following solutions can be expressed in terms of special functions.  Numerical libraries implement these special functions in a variety of ways.  In order to avoid confusion, we define the special functions here and indicate the syntax to produce them in a few libraries.

\begin{eqnarray}
\mathcal{G}(z) &=& \int_0^{\infty} t^{z-1} e^{-t} dt, 
\quad \text{complete gamma function}
\\
\mathcal{G}_{\rm L}(z,s) &=& \int_0^{s} t^{z-1} e^{-t} dt, 
\quad \text{lower incomplete gamma function}
\\
\mathcal{G}_{\rm U}(z,s) &=& \int_s^{\infty} t^{z-1} e^{-t} dt, 
\quad \text{upper incomplete gamma function}
\\
E_n(z) &=& \int_1^{\infty} t^{-n} e^{-z\,t} dt, 
\quad \text{exponential integral function}
\\
\mathcal{G}_{\rm L}(z,s) + \mathcal{G}_{\rm U}(z,s) &=& \mathcal{G}(z)
\\
E_n(z) &=& z^{n-1} \mathcal{G}_{\rm U}(1-n,z)
\end{eqnarray}
%
The complete gamma function $\mathcal{G}$ is usually designated $\Gamma$, but we reserve that symbol for photoionization rate.  Almost all libraries implement it as above.  However, it will not always be clear if a particular library has implemented the upper or lower version of the incomplete gamma function.  In addition, sometimes the incomplete gamma function is normalized by the complete gamma function.  Below we list the details of a few implementations.  


\subsubsection{Python - \texttt{mpmath}}
In the python module \texttt{mpmath} we can implement all of the functions above.  Assuming the line \texttt{import mpmath as mpm} is somewhere in the python script,
\begin{eqnarray}
\mathcal{G}(z) &=&  \texttt{mpm.gamma(z)}
\\
\mathcal{G}_{\rm L}(z,s) &=&  \texttt{mpm.gammainc(z,0,s)}
\\
\mathcal{G}_{\rm U}(z,s) &=&  \texttt{mpm.gammainc(z,s)}
\\
E_n(z) &=& \texttt{mpm.expint(n,z)}
\end{eqnarray}  


\subsubsection{Python - \texttt{scipy.special}}
In the python module \texttt{mpmath} we can implement all of the functions above.  Assuming the line \texttt{import scipy.special as ssp} is somewhere in the python script,
\begin{eqnarray}
\mathcal{G}(z) &=&  \texttt{ssp.gamma(z)}
\\
\mathcal{G}_{\rm L}(z,s) &=&  
\texttt{ssp.gammainc(z,s)} \times \texttt{ssp.gamma(z)}
\\
\mathcal{G}_{\rm U}(z,s) &=&  
\texttt{ssp.gammaincc(z,s)} \times \texttt{ssp.gamma(z)}
\\
E_n(z) &=& z^{n-1} \times
\texttt{ssp.gammaincc(1-n,z)} \times \texttt{ssp.gamma(1-n)}
\end{eqnarray}  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Basic Equations}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The evolution of the hydrogen neutral fraction, $\xHI = \nHI/\nH$, or equivalently ionized fraction, $\xHII = \nHII/\nH$, is a balance between photo/collisional ionizations and recombinations.

\begin{eqnarray}
\frac{d\xHI}{dt} &=& -(\Gamma + \gamma \, \nel) \xHI + 
\alpha \, \nel (1 - \xHI) \\
\frac{d\xHII}{dt} &=& (\Gamma + \gamma \, \nel)(1 - \xHII) - 
\alpha \, \nel \, \xHII
\end{eqnarray}
where $\Gamma$ is the photoionization rate, $\gamma$ is the collisional ionization rate, $\alpha$ is the recombination rate, and $\nel$ is the free electron number density.  If we decompose the free electron number density into a contribution from ionized hydrogen and a contribution from all heavier elements, we can write $\nel = (\xHII + y) \nH$.
The photoionization rate is in general given by, 
\begin{eqnarray}
\Gamma = \int d\Omega \int_{\nuth}^{\numax}
\frac{\Inu \signu}{h \nu} e^{-\taunu(\Omega)} d\nu = 
\ln 10 \int d\Omega \int_{\log \nuth}^{\log \numax}
\frac{\Inu \signu}{h} e^{-\taunu(\Omega)} d\log\nu
\end{eqnarray}
where $\Inu$ is the specific intensity of the radiation field and has cgs units of ergs s$^{-1}$ Hz$^{-1}$ cm$^{-2}$ Sr$^{-1}$, $\signu$ is the photoionization cross-section, and $\taunu(\Omega)$ is the frequency dependent optical depth along a given direction. If we model the spectrum of the radiation field as a power law and approximate the cross-section as a power law as well we can write, 

\begin{eqnarray}
w = \frac{\nu}{\nuth}, \quad
\Inu = I_0 w^{-\alpha}, \quad
\signu = \sigma_0 w^{-3}
\end{eqnarray}
If we consider a uniform radiation field then the angular integral can be performed and we can write, 

\begin{equation}
\Gamma = 4 \pi \frac{I_0 \sigma_0}{h} \int_{1}^{q}
w^{-(4+\alpha)} e^{-\tau_0 w^{-3}} dw, \quad
\tau_0 = \NHI \sigma_0 
\end{equation}
%
Concentrating on the integral we define, 
%
\begin{equation}
\mathcal{I}(q,\alpha,\tau_0) =
\int_{1}^{q} w^{-(4+\alpha)} e^{-\tau_0 w^{-3}} dw 
\end{equation}
where $q$ is the cut off frequency in Rydbergs.  We can write the solution to this integral in terms of well known special functions.  

In the special case of $q = \infty$ we can write,
\begin{equation}
\mathcal{I}(\alpha,\tau_0) = 
\frac{\mathcal{G}(A) - \mathcal{G}(A,\tau_0)}{3\tau_0^{A}}, 
\quad A = 1 + \frac{\alpha}{3}
\end{equation}
%
Otherwise,
\begin{equation}
\mathcal{I}(q,\alpha,\tau_0) = 
\frac{q^{-(3+\alpha)} E_{-\alpha/3}(\tau_0/q^3) - E_{-\alpha/3}(\tau_0)  }{3}
\end{equation}



\section{Exact Solutions}

Eqs. 1 and 2 can both be put into the form of a Riccati equation, 
\begin{eqnarray}
\frac{dx}{dt} = R x^2 + Q x + P.
\end{eqnarray}
Below we will determine the coefficients $R, Q, $ and $P$ for the neutral and ionized fractions. 


\subsection{Neutral Fraction Coefficients}

Grouping the terms in powers of $\xHI$, the first term in Eq. 1 becomes,
\begin{eqnarray} 
-(\Gamma + \gamma \nel) \xHI = \\
-\Gamma \xHI - \gamma [(1-\xHI) + y] \nH \xHI = \\
-\Gamma \xHI -\gamma \nH \xHI + \gamma \nH \xHI^2 - \gamma \nH y \xHI = \\
\gamma \nH \xHI^2  
- (\Gamma + \gamma \nH + \gamma \nH y ) \xHI
\end{eqnarray}
and the second becomes,

\begin{eqnarray} 
\alpha \nel (1 - \xHI) = \\
\alpha [(1-\xHI) + y] \nH (1-\xHI) = \\
\alpha \nH (1-\xHI)^2 + \alpha \nH y (1-\xHI) = \\
\alpha \nH (1 - 2\xHI + \xHI^2) + \alpha \nH y - \alpha \nH y \xHI =\\
\alpha \nH \xHI^2 - (2 \alpha \nH + \alpha \nH y) \xHI + 
\alpha \nH (1 +y)
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
(\gamma + \alpha) \nH \xHI^2 - 
\left[ \Gamma + \gamma \nH + \gamma \nH y  
+ 2 \alpha \nH + \alpha \nH y \right]  \xHI + 
\alpha \nH (1+y) = 0 \\
(\gamma + \alpha) \nH \xHI^2 -
\left[ \Gamma + (\gamma + 2\alpha) \nH + 
( \gamma + \alpha ) \nH y 
\right] \xHI + 
\alpha \nH (1+y) = 0 
\end{eqnarray}
Matching coefficients we get, 

\begin{eqnarray}
R &=& (\gamma + \alpha) \nH  \\
Q &=& -\left[ \Gamma + (\gamma + 2\alpha) \nH + 
(\gamma + \alpha) \nH y \right]  \\
  &=& -\left[ \Gamma + \alpha \nH + R (1+y) \right]  \\
P &=& \alpha \nH (1+y)  
\end{eqnarray}
Note that in the case of the neutral fraction, all the coefficients have a definite sign no matter what value the parameters take.  To gain some insight into 


\subsection{Ionized Fraction Coefficients}


Grouping the terms in powers of $\xHII$, the first term in Eq. 2 becomes,

\begin{eqnarray} 
(\Gamma + \gamma \nel) (1-\xHII) = \\
\left[ \Gamma + \gamma (\xHII + y) \nH \right]  (1-\xHII)  = \\
\Gamma + \gamma \nH \xHI + \gamma \nH y  -
\Gamma \xHI - \gamma \nH \xHI^2 - \gamma \nH y \xHI = \\
- \gamma \nH \xHI^2 + 
( \gamma \nH - \gamma \nH y ) \xHI
+ \Gamma + \gamma \nH y  
\end{eqnarray}
and the second becomes

\begin{eqnarray} 
- \alpha \nel \xHII = \\
- \alpha (\xHII + y) \nH \xHII = \\
-\alpha \nH \xHII^2 - \alpha \nH y \xHII
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
-( \gamma \nH + \alpha \nH ) \xHII^2 + 
( \gamma \nH - \Gamma - \gamma \nH y - \alpha \nH y) \xHII
+ \Gamma + \gamma \nH y = \\
-( \gamma + \alpha ) \nH \xHII^2 + 
\left[ \gamma \nH - \Gamma - (\gamma + \alpha) \nH y \right] \xHII
+ \Gamma + \gamma \nH y 
\end{eqnarray} 
Matching coefficients we get, 

\begin{eqnarray}
  R &=& -(\gamma + \alpha) \, \nH \\
  Q &=& \gamma \nH - \Gamma - 
  (\gamma + \alpha) \, \nH \, y  \nonumber \\
    &=& \gamma \nH - \Gamma + R \, y  \\
  P &=& \Gamma + \gamma \, \nH \, y   
\end{eqnarray}
In this case, the sign of $Q$ depends on the choice of parameters. 



\subsection{Solutions}

Now we have defined the coefficients $R$, $Q$, and $P$, we will explore solutions for $\xHI$ and $\xHII$.  Note that the following applies to both solutions (neutral and ionized fraction).  The only thing that changes is the values $R$, $Q$, and $P$.   

\subsubsection{Approach 1}
We first note that any quadratic expression can be factored in terms of its roots $x_{+}$ and $x_{-}$
\begin{eqnarray}
\frac{dx}{dt} = R x^2 + Q x + P = R ( x - x_{+} ) ( x - x_{-} ) 
\end{eqnarray}
where, 
\begin{eqnarray}
x_{\pm} = \frac{-Q \pm (Q^2 - 4 R P)^{1/2}}{2R}, \quad {\rm and} \quad
\delta x_\pm \equiv x_{+} - x_{-} = 
\frac{ (Q^2 - 4 R P)^{1/2}}{R}
\end{eqnarray}
The physical equilibrium solution is $x_{\rm eq} = x_{-}$.  To find the time dependent solution we integrate.

\begin{eqnarray}
\int_{x_0}^x \frac{dx'}{R ( x' - x_{+} ) ( x' - x_{-} ) } &=& \int_0^t dt 
\end{eqnarray}

\subsubsection{Method 1 - Partial Fractions}

\begin{eqnarray}
\frac{1}{(x-x_{+})(x-x_{\rm eq})} &=& 
\frac{A}{x-x_{+}} + \frac{B}{x-x_{\rm eq}} 
\\
1 &=&
A (x-x_{\rm eq}) + B (x-x_{+})
\end{eqnarray}
These relations must hold for any value of $x$.

\begin{eqnarray}
{\rm for} \quad x = x_{+}, &&\quad A = \frac{1}{x_{+}-x_{\rm eq}} 
= \frac{1}{\delta x_\pm} 
\\
{\rm for} \quad x = x_{\rm eq}, &&\quad B = \frac{1}{x_{\rm eq}-x_{+}} = -A
\end{eqnarray}
therefore, 

\begin{eqnarray}
  \frac{1}{(x-x_{+})(x-x_{\rm eq})} = 
  A \left( \frac{1}{x-x_{+}} - \frac{1}{x-x_{\rm eq}} \right)
\end{eqnarray}

\begin{eqnarray}
  \int_{x_0}^x \frac{dx'}{R ( x' - x_{+} ) ( x' - x_{-} ) } &=& \int_0^t dt 
  \\
  \frac{A}{R} \left( 
    \int_{x_0}^x \frac{dx'}{x' - x_{+}} - 
    \int_{x_0}^x \frac{dx'}{x' - x_{\rm eq}} 
  \right)  &=& t
  \\
  \left. \ln(x-x_{+}) \right|_{x_0}^x -  
  \left. \ln(x-x_{\rm eq}) \right|_{x_0}^x -  
   &=& \frac{R t}{A}
   \\
   \ln \left( \frac{ x-x_{+} }{ x_0 - x_{+} } \right) - 
   \ln \left( \frac{ x-x_{\rm eq} }{ x_0 - x_{\rm eq} } \right) 
   &=& R \delta x_{\pm} t
   \\
   \ln \left[ \left(  \frac{ x-x_{+} }{ x_0 - x_{+} } \right) 
   \left( \frac{ x_0 - x_{\rm eq} }{ x-x_{\rm eq} } \right) \right] 
   &=& R \delta x_{\pm} t
   \\
   \left(  \frac{ x-x_{+} }{ x_0 - x_{+} } \right) 
   \left( \frac{ x_0 - x_{\rm eq} }{ x-x_{\rm eq} } \right) 
   &=& \exp( R \delta x_{\pm} t )
   \\
   \left( \frac{x - x_{\rm eq}}{x - x_{+}} \right)    
   \left( \frac{x_0 - x_{+}}{x_0 - x_{\rm eq}}  \right)
   &=& \exp \left( -R \delta x_{\pm} t \right)
   \\
   \left( x - x_{\rm eq} \right) \left(x_0 - x_{+} \right)      
   &=& \left( x - x_{+} \right) \left( x_0 - x_{\rm eq}\right) 
   \exp \left( -R \delta x_{\pm} t \right)
\end{eqnarray}

Defining some helpful notation,
\begin{eqnarray}
  t_{\rm eq} &=& (R \delta x_{\pm})^{-1} = (Q^2 - 4 R P)^{-1/2} \\
  \mathcal{F} &=& \exp \left( - t/t_{\rm eq} \right)
\end{eqnarray}

\begin{eqnarray}
  \left( x - x_{\rm eq} \right) \left(x_0 - x_{+} \right)      
  &=& \left( x - x_{+} \right) \left( x_0 - x_{\rm eq}\right) 
  \mathcal{F}
  \\
  x        \left(x_0 - x_{+} \right) - 
  x_{\rm eq} \left(x_0 - x_{+} \right)      
  &=& x    \left( x_0 - x_{\rm eq}\right) \mathcal{F} - 
  x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} 
  \\
  x \left[ \left(x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq}\right) \mathcal{F} \right]   
  &=& x_{\rm eq} \left(x_0 - x_{+} \right) -       
  x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} 
  \\
  x \left[ \left(x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq}\right) \mathcal{F} \right]   
  &=& x_{\rm eq} \left[ \left( x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq} \right) \mathcal{F} \right] + 
  \\
  && 
  x_{\rm eq} \left( x_0 - x_{\rm eq} \right) \mathcal{F} -       
  x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} 
  \\
  x &=& x_{\rm eq} +  
  \frac{ 
    x_{\rm eq} \left( x_0 - x_{\rm eq} \right) \mathcal{F} -       
    x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} }
  {
    \left(x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq}\right) \mathcal{F} }
  \\
  x &=& x_{\rm eq} +  
  \frac{ 
    \left( x_{\rm eq} - x_{+} \right) 
    \left( x_0 - x_{\rm eq} \right) \mathcal{F}  }
  {
   x_0 (1 - \mathcal{F} ) + \mathcal{F} x_{\rm eq} - x_{+} }
  \\
  x &=& x_{\rm eq} +  \left( x_0 - x_{\rm eq} \right)
  \frac{ 
    \left( x_{\rm eq} - x_{+} \right) 
     \mathcal{F}  }
  {
   x_0 (1 - \mathcal{F} ) + \mathcal{F} x_{\rm eq} - x_{+} }
\end{eqnarray}

Another approach is the following, 


\begin{eqnarray}
  \int_{x_0}^x \frac{dx'}{R ( x' - x_{+} ) ( x' - x_{-} ) } &=& \int_0^t dt 
  \\
  \left.  
    \frac{\log(x'- x_{+}) - \log(x'-x_{-}) }{\delta x_{\pm}}
  \right|_{x_0}^x 
  &=& R t 
  \\ 
  \left. \log \left( \frac{x'-x_{+}}{x'- x_{-}}  \right) \right|_{x_0}^x 
  &=& R \delta x_{\pm} t 
  \\ 
  \log \left( \frac{x   - x_{+}}{x   - x_{-}}  \right) -
  \log \left( \frac{x_0 - x_{+}}{x_0 - x_{-}}  \right)
  &=& R \delta x_{\pm} t 
  \\ 
  \log \left[ 
    \left( \frac{x   - x_{+}}{x   - x_{-}} \right)    
    \left( \frac{x_0 - x_{-}}{x_0 - x_{+}}  \right)
  \right]
  &=& R \delta x_{\pm} t 
  \\ 
  \left( \frac{x   - x_{+}}{x   - x_{-}} \right)    
  \left( \frac{x_0 - x_{-}}{x_0 - x_{+}}  \right)
  &=& \exp \left( R \delta x_{\pm} t \right)
\end{eqnarray}
Taking the reciprocal and rearranging we get: 
\begin{eqnarray}
\left( \frac{x   - x_{-}}{x   - x_{+}} \right)    
\left( \frac{x_0 - x_{+}}{x_0 - x_{-}}  \right)
 &=& \exp \left( -R \delta x_{\pm} t \right)
\\ 
\left( \frac{x   - x_{-}}{x   - x_{+}} \right)    
 &=& \left( \frac{x_0 - x_{-}}{x_0 - x_{+}}  \right)
\exp \left( -R \delta x_{\pm} t \right) 
\end{eqnarray} 
Replacing $x_{-}$ with $x_{\rm eq}$ and defining some helpful notation,
\begin{eqnarray}
\mathcal{D} &=& \frac{x_0 - x_{\rm eq}}{x_0 - x_{+}} \\
t_{\rm eq} &=& (R \delta x_{\pm})^{-1} = (Q^2 - 4 R P)^{-1/2} \\
\mathcal{F} &=& \exp \left( - t/t_{\rm eq} \right)
\end{eqnarray}
we can write

\begin{eqnarray}
\left( \frac{x - x_{\rm eq}}{x - x_{+}} \right) &=& \mathcal{D} \mathcal{F}
\\
x - x_{\rm eq} &=& \left( x - x_{+} \right) \mathcal{D} \mathcal{F}
\\ 
x - x \mathcal{D} \mathcal{F} &=& 
x_{\rm eq} - x_{+} \mathcal{D} \mathcal{F}
\\ 
x \left( 1 - \mathcal{D} \mathcal{F} \right) &=& 
x_{\rm eq} - x_{+} \mathcal{D} \mathcal{F}
\\ 
x(t) &=& 
\frac{ x_{\rm eq} - x_{+} \mathcal{D} \mathcal{F}(t) }
{1 - \mathcal{D} \mathcal{F}(t) }
\end{eqnarray} 
This form makes it clear that as $t \rightarrow \infty$ and $\mathcal{F}(t) \rightarrow 0$, the solution $x(t) \rightarrow x_{\rm eq}$.  If we expand $\mathcal{D}$ we can present a form that makes it clear that as $t \rightarrow 0$ and $\mathcal{F}(t) \rightarrow 1$, the solution $x(t) \rightarrow x_0$.
\begin{eqnarray}
x &=& 
\frac{ x_{\rm eq}  - x_{+} \frac{x_0 - x_{\rm eq}}{x_0 - x_{+}} \mathcal{F} }
{1 - \frac{x_0 - x_{\rm eq}}{x_0 - x_{+}} \mathcal{F} }
\\
x &=& 
\frac{ x_{\rm eq} (x_0 - x_{+}) - x_{+} (x_0 - x_{\rm eq}) \mathcal{F} }
{(x_0 - x_{+}) - (x_0 - x_{\rm eq})  \mathcal{F} }
\\
x &=& 
\frac{ 
x_{\rm eq} x_0 - x_{\rm eq} x_{+} - x_{+} x_0 \mathcal{F} + x_{+} x_{\rm eq} \mathcal{F} }
{x_0 - x_{+} - x_0 \mathcal{F} + x_{\rm eq} \mathcal{F} }
\\
x(t) &=& 
\frac{ 
x_0 \left[ x_{\rm eq}  - x_{+} \mathcal{F}(t) \right]
- x_{\rm eq} x_{+} [1-\mathcal{F}(t)  ] }
{ [ x_{\rm eq} \mathcal{F}(t) - x_{+}] + x_0 [1- \mathcal{F}(t) ]  }
\end{eqnarray}

\subsubsection{Approach 2}

The time dependent solution of this Ricatti equation is, 
\begin{eqnarray}
 x(t) = \frac{1}{2R} \left[    
   -Q - d \frac{\tanh(t/t_{\rm c}) - x_{\rm s}}
   {1 - x_{\rm s} \tanh(t/t_{\rm c})}
   \right] \nonumber \\
  t_{\rm c} = \frac{2}{d}, \quad 
  x_{\rm s} = \frac{Q + 2 R \xo}{d}, \quad
  \xo = x(0), \quad
  d = \sqrt{Q^2 - 4 R P}
\end{eqnarray}
and the physical equilibrium solution is given by the negative root of
the quadratic formula.  

\begin{equation}
x^{\rm eq} = \frac{-Q -\sqrt{Q^2 - 4 R P}}{2R}
\end{equation}
In a gas composed only of Hydrogen and Helium, $y$ is bound between
zero and (1-$X$)/(2$X$) where $X$ is the Hydrogen mass fraction of the
gas.  For this solution then the ionization state is a function of
four variables $\nH$, $T$, $\Gamma$, and $y$.  This solution is best
evaluated numerically using the following form, 
\begin{equation}
x^{\rm eq} = \frac{P}{q}, \quad {\rm with} \quad
q = -\frac{1}{2} \left[ Q + {\rm sgn}(Q) \sqrt{Q^2 - 4RP} \right]
\end{equation}

\subsubsection{Physical Intuition}

We can gain some physical intuition about the parameters by considering some simplified cases.  Examining the neutral fraction has the great benefit that the coefficients $R$, $Q$, and $P$ have fixed signs.  In a pure hydrogen gas ($y=0$) with with no incident radiation ($\Gamma=0$) we have (for the neutral fraction case), 

\begin{eqnarray}
P = \alpha \nH  > 0, \quad 
R = P + \gamma \nH > P, \quad 
-Q = R + P
\end{eqnarray}
which makes it clear that $P < R < -Q$.  It follows that 
\begin{eqnarray}
Q^2 &=& (R+P)^2 \\
Q^2 &=& R^2 + P^2 + 2RP \\
Q^2 - 4RP &=& R^2 + P^2 - 2RP \\
Q^2 - 4RP &=& (R - P)^2
\end{eqnarray}
As $(R-P)^2$ must be positive, this implies that $Q^2 > 4RP$ and that $x_{+} > x_{-}$.  If we consider the case with a finite photoionization rate $\Gamma$ then the only coefficient that changes is $Q \rightarrow Q + \Gamma$.  If $Q^2 > 4RP$ is true then so is $(Q+\Gamma)^2 > 4RP$. 

\section{Numerical Solutions}

In the above solutions, $\Gamma$ refers to the photoionization rate 
at a single point.  In numerical work, discretization is a
necessity.  For infinitesimal fluid elements there is negligible
variation of the photo-ionization rate over the element and the above
equations are valid.  The larger the discretized fluid elements
become, the larger the errors associated with using a single $\Gamma$
become. To account for this, it is useful to consider seperately the 
contributions to the optical depth from the current fluid element 
(whose ionization state is being solved for) and all other fluid
elements.  In this case we replace $\Gamma$ with $\Gamma e^{-\tau_0
  \xHI}$ where $\Gamma$ is the photoionization rate calculated from
the optical depth along the ray \emph{excluding} the current fluid element
and $\tau_0$ is the current particle's contribution to the optical
depth assuming it is fully neutral.  There are three limiting cases we
want to examine, 
$\tau_0 \xHI << 1$, $\tau_0 \xHI >> 1$, and
$\tau_0 \xHI \approx \tau_0 \xHII \approx 1$.   


\subsubsection{Coefficients for $\tau_0 \xHI << 1$}

In this case we can expand the exponential term as 
$\Gamma e^{-\tau_0 \xHI} \approx \Gamma (1 - \tau_0 \xHI)$.  The first
term in Eq. 1 then becomes

\begin{eqnarray} 
-(\Gamma + \gamma \nel) \xHI = \\
-\Gamma (1-\tau_0 \xHI) \xHI - \gamma [(1-\xHI) + y] \nH \xHI = \\
-\Gamma \xHI + \Gamma \tau_0 \xHI^2  
-\gamma \nH \xHI + \gamma \nH \xHI^2 - \gamma \nH y \xHI = \\
(\Gamma \tau_0 + \gamma \nH) \xHI^2  
- (\Gamma + \gamma \nH + \gamma \nH y ) \xHI
\end{eqnarray}
and the second becomes

\begin{eqnarray} 
\alpha \nel (1 - \xHI) = \\
\alpha [(1-\xHI) + y] \nH (1-\xHI) = \\
\alpha \nH (1-\xHI)^2 + \alpha \nH y (1-\xHI) = \\
\alpha \nH (1 - 2\xHI + \xHI^2) + \alpha \nH y - \alpha \nH y \xHI =\\
\alpha \nH \xHI^2 - (2 \alpha \nH + \alpha \nH y) \xHI + 
\alpha \nH (1 +y)
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
\left[ \Gamma \tau_0 + (\gamma + \alpha) \nH \right] \xHI^2 - 
\left[ \Gamma + \gamma \nH + \gamma \nH y  
+ 2 \alpha \nH + \alpha \nH y \right]  \xHI + 
\alpha \nH (1+y) = 0 \\
\left[ \Gamma \tau_0 +(\gamma + \alpha) \nH \right] \xHI^2 -
\left[ \Gamma + (\gamma + 2\alpha) \nH + 
( \gamma + \alpha ) \nH y 
\right] \xHI + 
\alpha \nH (1+y) = 0 
\end{eqnarray}
\begin{eqnarray}
R &=& \Gamma \tau_0 +(\gamma + \alpha) \nH \\
Q &=& -\left[ \Gamma + (\gamma + 2\alpha) \nH + 
( \gamma + \alpha ) \nH y \right] \\
P &=& \alpha \nH (1+y)
\end{eqnarray}
and the physical solution is the negative root (c/q from Numerical
Recipes) of the quadratic formula.


\subsubsection{Coefficients for $\tau_0 \xHI >> 1$}

In this case we can hope that $\tau_0 (1-\xHI) = \tau_0 \xHII << 1$
and use the same procedure as above in Eq. 2.  This isn't guaranteed
to be true, so one should explicitly check the value of $\tau_0 \xHII$
when using this solution.  In this case, the exponential term is 
$\Gamma e^{-\tau_0 (1-\xHII)} = 
\Gamma e^{-\tau_0} e^{\tau_0 \xHII} \approx 
\Gamma e^{-\tau_0} (1 + \tau_0 \xHII)$.  The first
term in Eq. 2 then becomes

\begin{eqnarray} 
(\Gamma + \gamma \nel) (1-\xHII) = \\
\left[ \Gamma e^{-\tau_0} (1+\tau_0 \xHII) + 
\gamma (\xHII + y) \nH \right]  (1-\xHII)  = \\
\Gamma e^{-\tau_0} + \Gamma \tau_0 e^{-\tau_0} \xHII 
 + \gamma \nH \xHII + \gamma \nH y  -
\Gamma e^{-\tau_0} \xHII - \Gamma \tau_0 e^{-\tau_0} \xHII^2 
 - \gamma \nH \xHII^2 - \gamma \nH y \xHII = \\
-( \Gamma \tau_0 e^{-\tau_0} + \gamma \nH )\xHII^2 + 
( \Gamma \tau_0 e^{-\tau_0} - \Gamma e^{-\tau_0} 
+\gamma \nH - \gamma \nH y) \xHII
+ \Gamma e^{-\tau_0} + \gamma \nH y 
\end{eqnarray}
and the second becomes

\begin{eqnarray} 
- \alpha \nel \xHII = \\
- \alpha (\xHII + y) \nH \xHII = \\
-\alpha \nH \xHII^2 - \alpha \nH y \xHII
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
-( \Gamma \tau_0 e^{-\tau_0} + \gamma \nH + \alpha \nH ) \xHII^2 + 
( \Gamma \tau_0 e^{-\tau_0} - \Gamma e^{-\tau_0} 
+ \gamma \nH - \gamma \nH y 
- \alpha \nH y) \xHII
+ \Gamma e^{-\tau_0} + \gamma \nH y = \\
-\left[ \Gamma \tau_0 e^{-\tau_0} + (\gamma + \alpha) \nH \right] \xHII^2 + 
\left[ \Gamma e^{-\tau_0} (\tau_0 - 1) + \gamma \nH 
- (\gamma + \alpha) \nH y \right] \xHII 
+ \Gamma e^{-\tau_0} + \gamma \nH y = 0 
\end{eqnarray}
\begin{eqnarray}
R &=& -\left[ \Gamma \tau_0 e^{-\tau_0} + (\gamma + \alpha) \nH
  \right] \\
Q &=& \Gamma e^{-\tau_0} (\tau_0 - 1) 
+ \gamma \nH - (\gamma + \alpha) \nH y \\
P &=& \Gamma e^{-\tau_0} + \gamma \nH y
\end{eqnarray} 


\subsubsection{Solution for $\tau_0 \xHI \approx \tau_0 \xHII \approx 1$}

In this case we cannot expand the exponential term and we do not get a
quadratic equation for $\xHI$.  The first term in Eq. 1 then becomes,

\begin{eqnarray} 
-(\Gamma + \gamma \nel) \xHI = \\
-\Gamma e^{-\tau_0 \xHI} \xHI - \gamma [(1-\xHI) + y] \nH  \xHI = \\
-\Gamma e^{-\tau_0 \xHI} \xHI - \gamma \nH \xHI + \gamma \nH \xHI^2  
- \gamma \nH y \xHI = \\
\gamma \nH \xHI^2 - ( \gamma \nH + \gamma \nH y) \xHI 
- \Gamma e^{-\tau_0 \xHI} \xHI
\end{eqnarray}
and the second becomes

\begin{eqnarray} 
\alpha \nel (1 - \xHI) = \\
\alpha [(1-\xHI) + y] \nH (1-\xHI) = \\
\alpha \nH (1-\xHI)^2 + \alpha \nH y (1-\xHI) = \\
\alpha \nH (1 - 2\xHI + \xHI^2) + \alpha \nH y - \alpha \nH y \xHI =\\
\alpha \nH \xHI^2 - (2 \alpha \nH + \alpha \nH y) \xHI + 
\alpha \nH (1 +y)
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
(\gamma + \alpha) \nH \xHI^2 - 
\left[ (\gamma + 2\alpha) \nH + (\gamma + \alpha) \nH y \right] \xHI + 
\alpha \nH (1+y) - \Gamma e^{-\tau_0 \xHI} \xHI 
\end{eqnarray}
\begin{eqnarray}
\frac{d\xHI}{dt} &=& R \xHI^2 + Q \xHI + P + S  
e^{-\tau_0 \xHI} \xHI  = f(\xHI) = 0\nonumber \\
 R  &=& (\gamma + \alpha) \nH \\
 Q  &=& -\left[ (\gamma + 2\alpha) \nH + (\gamma + \alpha) \nH y
   \right] \\
 P  &=&  \alpha \nH (1+y) \\ 
 S  &=& -\Gamma
\end{eqnarray}
for an iterative solution it is useful to have the derivative,

\begin{eqnarray}
\frac{df}{d\xHI} = 2 R \xHI + Q + S e^{-\tau_0 \xHI}
(1 - \xHI \tau_0)
\end{eqnarray}

For certain cases, care must be taken when evaluating the
terms numerically.  The exponential term $e^{-\tau_0 \xHI}$ can easily
cause underflows.  It is also important to keep in mind that with
non-zero recombination and collisional ionization rates, the extreme
values $\xHI = 1$, and $\xHI = 0$ do not occur and will not solve the
above equation.  In addition, we note that this solution is equivalent
to the analytic neutral fraction solution in the case of zero
photo-ionization rate.  This equation can have multiple physical
solutions for some sets of parameters.  The approach that makes the
most sense to me is to begin with the physical solution in the zero
photo-ionization rate case and slowly increase $\Gamma$, tracking the
physical solution. 



\begin{figure}
  \begin{center}
    \includegraphics[width=0.45\textwidth]{images/atomic_rates}
  \end{center}
  \caption{Recombination and collisional ionization rates as a
  function of temperature.  Lower panel: The red (blue) dashed (solid)
  line in the bottom panel is the case B (case A) recombination rate
  while the green dash-dot line is the collisional ionization rate 
  \cite{1997MNRAS.292...27H}.  Top panel: The ratio of the collisional
  ionization rate to both recombination rates. }
\end{figure}





\section{Plane Parallel Radiation on Uniform Slab w/o Collisions}

Our first problem is that of monochromatic plane parallel radiation 
incident on a semi-infinite slab of pure Hydrogen ($y = 0$)with an 
arbitrary density vs. depth relationship $\nH = \nH(z)$.  
If we choose the temperature of this gas to be low enough, we can also
set $\gamma$ to zero.  We choose the coordinate system such that
the surface of the slab is coincident with the $x-y$ plane and the
positive $z$ direction extends into the slab. Our goal is to calculate
the ionization fraction $\xHII(z)$ as a function of $z$ given a flux
of incoming photons per second per cm$^2$, $F_0$.  

The photoionization rate is related to 
$\tau(z) = \int_0^z \nH \xHI \sigth ds$ as, 
$\Gamma(z) = F(z) \sigth = F_0 \sigth e^{-\tau(z)} $. At the
surface of the slab, we have the boundary condition, 
$\tau(0) = 0$, and so $\Gamma(0) = \Gamma_0 = F_0 \sigth$.
With $\gamma$ and $y$ equal to zero, the equilibrium neutral/ionization
fractions obey the equations,
\begin{eqnarray}
\frac{\xHI}{(1-\xHI)^2} = 
\frac{\alpha \nH}{\Gamma} = 
\frac{\alpha \nH}
{\Gamma_0 {\rm Exp}[{-\int_0^z \nH \xHI \sigth ds}]}\\
\frac{\xHII^2}{1-\xHII} = \frac{\Gamma}{\alpha \nH} = 
\frac{\Gamma_0 {\rm Exp}[{-\int_0^z \nH (1-\xHII) \sigth ds}]}{\alpha \nH}
\end{eqnarray}
The values at $z=0$ are,

\begin{eqnarray}
\frac{\xHIzero}{(1-\xHIzero)^2} = 
\frac{\alpha \nHzero}{\Gamma_0} &=& \frac{1}{\mathcal{C}} > 0, \quad
\xHIzero^2 - (2+\mathcal{C}) \xHIzero + \mathcal{C} = 0\\
\frac{\xHIIzero^2}{1-\xHIIzero} = 
\frac{\Gamma_0}{\alpha \nHzero} &=& \mathcal{C} > 0, \quad
\xHIIzero^2 + \mathcal{C} \xHIIzero - \mathcal{C} = 0
\end{eqnarray}
The physical roots of the quadratic formula are negative (positive)
for the neutral (ionized) fraction. 

\begin{eqnarray}
\xHIzero = \frac{2+\mathcal{C} - \sqrt{\mathcal{C}^2 + 4\mathcal{C}}}{2}\\
\xHIIzero = \frac{-\mathcal{C} + \sqrt{\mathcal{C}^2 + 4\mathcal{C}}}{2}
\end{eqnarray}

\subsection{Neutral Fraction Details}

Lets start with solving for the $z$ dependent neutral 
fraction $\xHI$.  We have,

\begin{eqnarray}
\frac{\xHI}{(1-\xHI)^2} \mathcal{C} &=& 
{\rm Exp} \left[ \int_0^z \nH \xHI \sigth ds \right] 
\\
\ln \left[ \frac{\xHI}{(1-\xHI)^2} \mathcal{C} \right] &=& 
\int_0^z \nH \xHI \sigth ds 
\\
\frac{d}{dz} \left[ \ln(\xHI) - 2\ln(1-\xHI) + \ln(\mathcal{C}) \right]&=& 
\frac{d}{dz} \int_0^z \nH \xHI \sigth ds 
\\
\frac{d}{dz} \left[ \ln(\xHI) - 2\ln(1-\xHI)  \right]&=& 
\nH \xHI \sigth 
\\
\frac{1}{\xHI} \frac{d\xHI}{dz} + \frac{2}{1-\xHI} \frac{d\xHI}{dz} &=&
\nH \xHI \sigth 
\\
\left[ \frac{1}{\xHI^2} + \frac{2}{\xHI(1-\xHI)} \right]
\frac{d\xHI}{dz} &=& \nH \sigth 
\\
\int_{\xHIzero}^x 
\left[ \frac{1}{\xHI^2} + \frac{2}{\xHI(1-\xHI)} \right]
d\xHI &=& 
\int_0^z \nH \sigth dz = \tau_{\rm H}
\\ 
-\frac{1}{\xHI} \Big|_{\xHIzero}^x +
2 \left[ \ln \xHI - \ln (\xHI-1) \right] \Big|_{\xHIzero}^x + 
&=& \tau_{\rm H}
\\
\left[ \frac{1}{\xHIzero} - \frac{1}{x}  \right] + 
2 \left[ \ln \frac{x}{\xHIzero} + \ln \frac{\xHIIzero-1}{x-1} \right]  
&=& \tau_{\rm H} 
\\
\left[ \frac{1}{\xHIzero} - \frac{1}{x}  \right] + 
2 \left[ \ln \frac{x(\xHIzero-1)}{\xHIzero(x-1)} \right]  
&=& \tau_{\rm H} 
\\
\left[ \frac{1}{\xHIzero} - \frac{1}{x}  \right] + 
2 \left[ \ln \frac{x(1-\xHIzero)}{\xHIzero(1-x)} \right]  
&=& \tau_{\rm H} 
\end{eqnarray}
and since $x$ is arbitrary we can write,

\begin{eqnarray}
\tau_{\rm H} = \left[ \frac{1}{\xHIzero} - \frac{1}{\xHI}  \right] + 
2  \ln \left[ \frac{\xHI(1-\xHIzero)}{\xHIzero(1-\xHI)} \right]
\end{eqnarray}

This relation gives the optical depth from the surface of the slab in terms
of the neutral fraction $\xHI$.  This is the inverse of the
function we were looking for, but it can be used to map out the solution.
Starting from $\xHI=\xHIzero$ we can increase $\xHI$ in as many
increments as we need to cover a suitable range in $\tau_{\rm H}$.  In
this case, the solution applies to arbitrary density vs distance
$\nH(z)$ relationships.

\subsection{Ionized Fraction Details}

To be pedagogical and to check our result, we can do the same
calculation for the ionized fraction $\xHII$.  We have,


\begin{eqnarray}
\frac{\xHII^2}{1-\xHII} \frac{1}{\mathcal{C}} &=& 
{\rm Exp} \left[ -\int_0^z \nH (1-\xHII) \sigth ds \right] \\
\ln \left[ \frac{\xHII^2}{1-\xHII} \frac{1}{\mathcal{C}} \right] &=& 
-\int_0^z \nH (1-\xHII) \sigth ds \\
\frac{d}{dz} \left[ \ln(\xHII^2) - \ln(1-\xHII) - \ln(\mathcal{C}) \right]&=& 
-\frac{d}{dz} \int_0^z \nH (1-\xHII) \sigth ds \\
\frac{d}{dz} \left[ 2\ln(\xHII) - \ln(1-\xHII)  \right]&=& 
-\nH (1-\xHII) \sigth \\
\frac{2}{\xHII} \frac{d\xHII}{dz} + \frac{1}{1-\xHII} \frac{d\xHII}{dz} &=&
-\nH (1-\xHII) \sigth \\
\left[ \frac{2}{\xHII (1-\xHII)} + \frac{1}{(1-\xHII)^2} \right]
\frac{d\xHII}{dz} &=& -\nH \sigth \\
\int_{\xHIIzero}^x 
\left[ \frac{2}{\xHII (1-\xHII)} + \frac{1}{(1-\xHII)^2} \right]
d\xHII &=& -\int_0^z \nH \sigth dz = - \tau_{\rm H}\\ 
2 \ln \frac{\xHII}{\xHII -1} \Big|_{\xHIIzero}^x + 
\frac{1}{(1-\xHII)} \Big|_{\xHIIzero}^x 
&=& -\tau_{\rm H}\\
2 \left[ \ln \frac{x}{x-1} - \ln \frac{\xHIIzero}{\xHIIzero-1} \right] + 
\left[ \frac{1}{1-x} - \frac{1}{1-\xHIIzero} \right] &=& -\tau_{\rm H} \\
-2 \left[ \ln \frac{x}{x-1} - \ln \frac{\xHIIzero}{\xHIIzero-1} \right] - 
\left[ \frac{1}{1-x} - \frac{1}{1-\xHIIzero} \right] &=& \tau_{\rm H} \\
-2 \ln \left( \frac{x}{x-1} \frac{\xHIIzero-1}{\xHIIzero} \right) - 
\left( \frac{1}{1-x} - \frac{1}{1-\xHIIzero} \right) &=& \tau_{\rm H} \\
2 \ln \left( \frac{1-x}{x} \frac{\xHIIzero}{1-\xHIIzero} \right) - 
\left( \frac{1}{1-x} - \frac{1}{1-\xHIIzero} \right) &=& \tau_{\rm H}
\end{eqnarray}
and since $x$ is arbitrary we can write,

\begin{eqnarray}
\tau_{\rm H} = 
\left[ \frac{1}{1-\xHIIzero} - \frac{1}{1-\xHII}   \right] +
2 \ln \left[ \frac{\xHIIzero(1-\xHII)}{\xHII(1-\xHIIzero)} \right]  
\end{eqnarray}

This relation gives the optical depth from the surface of the slab in terms
of the ionized fraction $\xHII$.  This is the inverse of the
function we were looking for, but it can be used to map out the solution.
Starting from $\xHII=\xHIIzero$ we can decrease $\xHII$ in as many
increments as we need to cover a suitable range in $\tau_{\rm H}$.  In
this case, the solution applies to arbitrary density vs distance
$\nH(z)$ relationships.


\section{Plane Parallel Radiation on Uniform Slab with Collisions}

Here we examine the case of monochromatic plane parallel radiation 
incident on a semi-infinite slab of pure Hydrogen ($y = 0$) with a 
constant density  $\nH$.  
In this derivation, we include the effects of collisional
ionizations. 
We choose the coordinate system such that the surface of the slab is  
coincident with the $x-y$ plane and the positive $z$ direction extends
into the slab. Our goal is to calculate the ionization state 
$\xHI(z) = 1-\xHII(z)$ as a function of $z$ given a flux of incoming
photons per
second per cm$^2$, $F_0$.  
The photoionization rate is related to 
$\tau(z) = \int_0^z \nH \xHI \sigth ds$ as, 
$\Gamma(z) = F(z) \sigth = F_0 \sigth e^{-\tau(z)} $. At the
surface of the slab, we have the boundary condition, 
$\tau(0) = 0$, and so $\Gamma(0) = \Gamma_0 = F_0 \sigth$.
The ionization states at $z=0$ ($\xHIzero$ and $\xHIIzero$) can be
calculated by first calculating
$R$, $Q$, and $P$ and then using the solutions from \S 3 above with
$\nH, T, \Gamma_0, y$.

\subsection{Neutral Fraction Solution}

Lets begin with the $\xHI$ solution, 

\begin{eqnarray}
\frac{d\xHI}{dt} &=& -(\Gamma + \gamma \, \nel) \xHI + 
\alpha \, \nel (1 - \xHI)
\end{eqnarray}
therefore, 

\begin{eqnarray}
(\Gamma + \gamma \nel) \xHI &=&   
\alpha \nel (1-\xHI) \nonumber 
\\
\Gamma \xHI + \gamma  \nH (1 - \xHI) \xHI  &=&   
\alpha \, \nH (1 - \xHI)^2 \nonumber 
\end{eqnarray}
therefore, 

\begin{eqnarray}
\frac{\Gamma}{\Gamma_0}  &=&   
\frac{\alpha \nH (1 - \xHI)^2 - \gamma \nH (1 - \xHI) \xHI }
{\Gamma_0 \xHI}
 \nonumber 
\\
{\rm Exp} \left[ -\int_0^z \nH \xHI \sigth ds \right]  &=& 
\frac{\alpha \nH (1 - \xHI)^2 - \gamma \nH (1 - \xHI) \xHI }
{\Gamma_0 \xHI} \nonumber 
\\ 
-\int_0^z \nH \xHI \sigth ds   &=& 
\ln \left\{ \frac{\alpha \nH (1 - \xHI)^2 - \gamma \nH (1 - \xHI) \xHI }
{\Gamma_0 \xHI} \right\} \nonumber 
\\
-\int_0^z \nH \xHI \sigth ds   &=& 
\ln \left[ \alpha \nH (1 - \xHI)^2 - \gamma \nH (1 - \xHI) \xHI \right] - 
\ln \left[ \Gamma_0 \xHI \right] \nonumber 
\\
-\int_0^z \nH \xHI \sigth ds   &=& 
\ln \left[ \nH \right] +
\ln \left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] - 
\ln \left[ \Gamma_0 \xHI \right] \nonumber 
\\
\int_0^z \nH \xHI \sigth ds   &=& 
-\ln \left[ \nH \right] 
-\ln \left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
+\ln \left[ \Gamma_0 \xHI \right] 
\nonumber 
\\
\frac{d}{dz}\int_0^z \nH \xHI \sigth ds   &=& 
\frac{d}{dz} \left\{ 
-\ln \left[ \nH \right] 
-\ln \left[ \alpha (1-\xHI)^2 - \gamma (1-\xHI) \xHI \right]
+\ln \left[ \Gamma_0 \right]  
+\ln \left[ \xHI \right] 
\right\} \nonumber 
\\
\nH \xHI \sigth &=& 
-\frac{1}{\nH} \frac{d\nH}{dz} - 
\frac{d}{dz}  \left\{ 
\ln \left[ \alpha (1-2\xHI+\xHI^2) 
-\gamma \xHI + \gamma \xHI^2 \right]
-\ln \left[ \xHI \right] 
\right\} \nonumber 
\\
\nH \xHI \sigth &=& 
-
\left\{
\frac{ \alpha (-2+2\xHI) - 
\gamma + 2 \gamma \xHI}
{\alpha (1-\xHI)^2 - \gamma (1-\xHI) \xHI} 
-\frac{1}{\xHI} 
\right\} \frac{d\xHI}{dz}
 \nonumber 
\\
\nH \xHI \sigth &=& 
- 
\left\{
\frac{ 2 (\alpha + \gamma) \xHI - 
2 \alpha - \gamma }
{\alpha (1-\xHI)^2 - \gamma (1-\xHI) \xHI} 
-\frac{1}{\xHI} 
\right\} \frac{d\xHI}{dz}
 \nonumber 
\\
\int_0^z \nH \sigth dz  &=& 
\int_{\xHIzero}^x \left\{ 
\frac{1}{\xHI^2} - 
\frac{ 2 (\alpha + \gamma) \xHI - 
2 \alpha - \gamma }
{\alpha (1-\xHI)^2 \xHI - \gamma (1-\xHI) \xHI^2} 
\right\} d\xHI \nonumber 
\\
\tau_{\rm H} &=& 
\frac{-1}{\xHI} \Big|_{\xHIzero}^x
- \int_{\xHIzero}^x 
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{\alpha (1-\xHI)^2 \xHI - \gamma (1-\xHI) \xHI^2} 
 d\xHI \nonumber
\\
\tau_{\rm H} &=& 
\left( \frac{1}{\xHIzero} - \frac{1}{x} \right)
- \mathcal{I} \nonumber
\end{eqnarray}
Examining the integrand alone for a moment we have, 

\begin{eqnarray}
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{\alpha (1-\xHI)^2 \xHI - \gamma (1-\xHI) \xHI^2} = 
% \\
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{\alpha (1-2\xHI+\xHI^2) \xHI - \gamma \xHI^2 + \gamma \xHI^3} = 
\nonumber
\\
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{(\alpha+\gamma)\xHI^3 - (2\alpha + \gamma) \xHI^2 + \alpha \xHI} = 
%  \\
\frac{2 \xHI + \frac{-2 \alpha - \gamma}{\alpha + \gamma} }
{\xHI^3 - \frac{2\alpha + \gamma}{\alpha + \gamma} \xHI^2 + 
\frac{\alpha}{\alpha + \gamma} \xHI} = 
\nonumber
\\
\frac{2 \xHI - \frac{2 \alpha + \gamma}{\alpha + \gamma} }
{\xHI \left( \xHI^2 - \frac{2\alpha + \gamma}{\alpha + \gamma} \xHI + 
\frac{\alpha}{\alpha + \gamma} \right)} = 
% \\
\frac{2 \xHI - \left(
\frac{\alpha + \gamma}{\alpha + \gamma} + 
\frac{\alpha}{\alpha + \gamma} \right)}
{\xHI \left[ \xHI^2 - 
\left( \frac{\alpha + \gamma}{\alpha + \gamma} +
\frac{\alpha}{\alpha + \gamma} \right) 
 \xHI + 
\frac{\alpha}{\alpha + \gamma} \right]} = 
\nonumber
\\
\frac{2 \xHI - \left(
1 + \frac{\alpha}{\alpha + \gamma} \right)}
{\xHI \left[ \xHI^2 - 
\left( 1 + \frac{\alpha}{\alpha + \gamma} \right)
 \xHI + 
\frac{\alpha}{\alpha + \gamma} \right]} = 
% \\
\frac{2 \xHI - (1+\mathcal{A}) }
{\xHI \left[ \xHI^2 - (1+\mathcal{A}) \xHI + \mathcal{A} \right]} = 
\nonumber
\\
\frac{2 \xHI - (1+\mathcal{A}) }
{\xHI (\xHI - 1)(\xHI- \mathcal{A})}, 
\quad {\rm with} \quad 
\mathcal{A} = \frac{\alpha}{\alpha+\gamma}
 = \xHI^{\rm ci,eq} > \xHI
\end{eqnarray}
where $\mathcal{A}$ is the equilibrium neutral fraction in the absence
of photo-ionizations.  Now we can decompose this into partial fractions, 

\begin{eqnarray}
\frac{2 \xHI - (1+\mathcal{A}) }
{\xHI (\xHI - 1)(\xHI- \mathcal{A})} &=& 
\nonumber
\\
\left\{ 2 \xHI - (1+\mathcal{A}) \right\}
\left\{ 
\frac{\mathcal{A}^{-1}}{\xHI} + 
\frac{(1-\mathcal{A})^{-1}}{\xHI-1} + 
\frac{[\mathcal{A}(\mathcal{A}-1)]^{-1}}{\xHI-\mathcal{A}} 
\right\} &=& 
\nonumber
\\
\frac{2\mathcal{A}^{-1}}{1} -
\frac{\mathcal{A}^{-1} + 1}{\xHI} + 
\frac{2\xHI(1-\mathcal{A})^{-1}}{\xHI-1} - 
\frac{(1+\mathcal{A})/(1-\mathcal{A})}{\xHI-1} + 
\frac{2\xHI[\mathcal{A}(\mathcal{A}-1)]^{-1}}{\xHI-\mathcal{A}} - 
\frac{(1+\mathcal{A})[\mathcal{A}(\mathcal{A}-1)]^{-1}}{\xHI-\mathcal{A}} 
&=& \nonumber
\\
\frac{2}{\mathcal{A}} -
\frac{\mathcal{A}^{-1} + 1}{\xHI} - 
\frac{2}{1-\mathcal{A}} \frac{\xHI}{1-\xHI} + 
\frac{1+\mathcal{A}}{1-\mathcal{A}} \frac{1}{1-\xHI} - 
\frac{2}{\mathcal{A}(\mathcal{A}-1)} \frac{\xHI}{\mathcal{A}-\xHI} + 
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\frac{1}{\mathcal{A}-\xHI} 
&&
\end{eqnarray}
In the last step I've swapped negative signs on some of the terms so
that the individual integrations yield non-imaginary values.  Now, we can
integrate these term by term,

\begin{eqnarray}
\frac{2}{\mathcal{A}} \int_{\xHIzero}^x d\xHI  = 
\frac{2}{\mathcal{A}} (x - \xHIzero)
\\
-(\mathcal{A}^{-1} + 1)\int_{\xHIzero}^x \frac{d\xHI}{\xHI} = 
-(\mathcal{A}^{-1} + 1) \ln({x/\xHIzero})
\\
\frac{-2}{1-\mathcal{A}} 
\int_{\xHIzero}^x 
\frac{\xHI d\xHI}{1-\xHI}  = 
\frac{-2}{1-\mathcal{A}} 
[ -\xHI - \ln(1-\xHI) ] \Big|_{\xHIzero}^x = \nonumber
\\
\frac{-2}{1-\mathcal{A}} 
\left[ (\xHIzero-x) + \ln \left( \frac{1-\xHIzero}{1-x} \right) \right]
\\
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\int_{\xHIzero}^x \frac{d\xHI}{1-\xHI} = 
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\ln \left( \frac{1-\xHIzero}{1-x} \right) 
\\
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\int_{\xHIzero}^x \frac{\xHI d\xHI}{\mathcal{A}-\xHI} = 
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\left\{ 
\mathcal{A}(1-\ln \mathcal{A}^{-1}) - 
\left[ \xHI + \mathcal{A} \ln(\mathcal{A} - \xHI)
\right] \right\} 
\Big|_{\xHIzero}^x = \nonumber
\\
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\left[ (\xHIzero-x) + 
\mathcal{A} 
\ln \left( 
\frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} 
 \right) \right]
\\
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\int_{\xHIzero}^x \frac{d\xHI}{\mathcal{A}-\xHI} = 
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\ln \left( 
\frac{\mathcal{A} - \xHIzero}{\mathcal{A}-x}
\right)
\end{eqnarray}
The integral $\mathcal{I}$ is then,

\begin{eqnarray}
\mathcal{I} = 
\frac{2}{\mathcal{A}} (x - \xHIzero) - 
(\mathcal{A}^{-1} + 1) \ln({x/\xHIzero}) + 
\nonumber
\\
\frac{-2}{1-\mathcal{A}} 
\left[ (\xHIzero-x) + 
\ln \left( \frac{1-\xHIzero}{1-x} \right) \right] + \nonumber
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\ln \left( \frac{1-\xHIzero}{1-x} \right) +
\\
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\left[ (\xHIzero-x) + 
\mathcal{A} 
\ln \left( 
\frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} 
 \right) \right] + 
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\ln \left( 
\frac{\mathcal{A} - \xHIzero}{\mathcal{A}-x}
\right) 
\end{eqnarray}
Rearranging negative signs such that terms in $()$ are positive,

\begin{eqnarray}
\mathcal{I} = \frac{2}{\mathcal{A}} (x - \xHIzero) - 
(1 + \mathcal{A}^{-1}) \ln({x/\xHIzero}) + 
\nonumber
\\
\frac{2}{1-\mathcal{A}} 
\left[ (x - \xHIzero) - 
\ln \left( \frac{1-\xHIzero}{1-x} \right) \right] + \nonumber
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\ln \left( \frac{1-\xHIzero}{1-x} \right) +
\\
\frac{-2}{\mathcal{A}(1-\mathcal{A})}
\left[ (x-\xHIzero) - 
\mathcal{A} 
\ln \left( 
\frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} 
 \right) \right] - 
\frac{1+\mathcal{A}}{\mathcal{A}(1-\mathcal{A})}
\ln \left( 
\frac{\mathcal{A} - \xHIzero}{\mathcal{A}-x}
\right) = 
\end{eqnarray}
Simplifying we see,

\begin{eqnarray}
\left[ 
\frac{2}{\mathcal{A}} + 
\frac{2}{1-\mathcal{A}} -
\frac{2}{\mathcal{A}(1-\mathcal{A})} 
\right] (x-\xHIzero) = 0
\\
\left[ 
-\frac{2}{1-\mathcal{A}} + 
\frac{1+\mathcal{A}}{1-\mathcal{A}}
\right] 
\ln \left( \frac{1-\xHIzero}{1-x} \right) = 
%-----------------
- 
\ln \left( \frac{1-\xHIzero}{1-x} \right)
\\
\left[ 
\frac{2\mathcal{A}}{\mathcal{A}(1-\mathcal{A})} -
\frac{1+\mathcal{A}}{\mathcal{A}(1-\mathcal{A})} 
\right] 
\ln \left( \frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} \right) = 
%-----------------------------
\left[ 
-\frac{1}{\mathcal{A}} 
\right] 
\ln \left( \frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} \right) 
\end{eqnarray}
and the total inverse solution is, 

\begin{eqnarray}
\tau_{\rm H} &=& 
\left( \frac{1}{\xHIzero} - \frac{1}{x} \right) +
(1 + \mathcal{A}^{-1}) 
\ln \left( \frac{x}{\xHIzero} \right) +
\ln \left( \frac{1-\xHIzero}{1-x} \right) + 
\frac{1}{\mathcal{A}} 
\ln \left( \frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} \right) 
\\
&=&
\left( \frac{1}{\xHIzero} - \frac{1}{x} \right) +
\ln \left[ \frac{x(1-\xHIIzero)}{\xHIzero(1-x)} \right] + 
\frac{1}{\mathcal{A}}
\ln \left[ \frac{x(\mathcal{A}-\xHIIzero)}{\xHIzero(\mathcal{A}-x)} \right] 
\end{eqnarray}
Finally, because $x$ is arbitrary, we can write,

\begin{eqnarray}
\tau_{\rm H} = 
\left( \frac{1}{\xHIzero} - \frac{1}{\xHI} \right) +
\ln \left[ \frac{\xHI(1-\xHIIzero)}{\xHIzero(1-\xHI)} \right] + 
\frac{1}{\mathcal{A}}
\ln \left[ 
\frac{\xHI(\mathcal{A}-\xHIIzero)}
{\xHIzero(\mathcal{A}-\xHI)} 
\right] 
\end{eqnarray}




\subsection{Ionized Fraction Solution}

Now lets move on to solving for the $z$ dependent ionization
fraction $\xHII$.  We have,

\begin{eqnarray}
\frac{d\xHII}{dt} &=& (\Gamma + \gamma \, \nel)(1 - \xHII) - 
\alpha \, \nel \, \xHII = 0
\end{eqnarray}
therefore, 
\begin{eqnarray}
(\Gamma + \gamma \, \nel)(1 - \xHII) &=&   
\alpha \, \nel \, \xHII \nonumber 
\\
\Gamma (1 - \xHII) + \gamma \, \nH \xHII (1 - \xHII) &=&   
\alpha \, \nH \, \xHII^2 \nonumber 
\end{eqnarray}
therefore, 
\begin{eqnarray}
\frac{\Gamma}{\Gamma_0}  &=&   
\frac{\alpha \, \nH \, \xHII^2 - \gamma \nH \xHII (1 - \xHII)}
{\Gamma_0(1-\xHII)}
 \nonumber 
\\
{\rm Exp} \left[ -\int_0^z \nH (1-\xHII) \sigth ds \right]  &=& 
\frac{\alpha \nH \xHII^2 - \gamma \nH \xHII (1-\xHII)}
{\Gamma_0(1-\xHII)} \nonumber 
\\ 
-\int_0^z \nH (1-\xHII) \sigth ds   &=& 
\ln \left\{ \frac{\alpha \nH \xHII^2 - \gamma \nH \xHII (1-\xHII)}
{\Gamma_0(1-\xHII)} \right\} \nonumber 
\\
-\int_0^z \nH (1-\xHII) \sigth ds   &=& 
\ln \left[ \alpha \nH \xHII^2 - \gamma \nH \xHII (1-\xHII) \right] - 
\ln \left[ \Gamma_0(1-\xHII) \right] \nonumber 
\\
\int_0^z \nH (1-\xHII) \sigth ds   &=& 
\ln \left[ \Gamma_0(1-\xHII) \right] - 
\ln \left[ \alpha \nH \xHII^2 - \gamma \nH \xHII (1-\xHII) \right] 
\nonumber 
\\
\frac{d}{dz}\int_0^z \nH (1-\xHII) \sigth ds   &=& 
\frac{d}{dz} \left\{ \ln \left[ \Gamma_0(1-\xHII) \right] - 
\ln \left[ \alpha \nH \xHII^2 - \gamma \nH \xHII (1-\xHII) \right]
\right\} \nonumber 
\\
\nH (1-\xHII) \sigth &=& 
\frac{d}{dz} \left\{ \ln \left[ \Gamma_0 \right] +
\ln \left[ (1-\xHII) \right] - 
\ln \left[ \alpha \xHII^2 - \gamma \xHII (1-\xHII) \right] - 
\ln \left[ \nH \right]
\right\} \nonumber 
\\
\nH (1-\xHII) \sigth &=& 
\left\{ 
\frac{-1}{1-\xHII}   - 
\frac{2 \alpha \xHII - \gamma + 2 \gamma \xHII}
{\alpha \xHII^2 - \gamma \xHII (1-\xHII)} 
\right\} \frac{d\xHII}{dz} \nonumber 
\\
\int_0^z \nH \sigth dz &=& 
\int_{\xHIIzero}^x \left\{ 
\frac{-1}{(1-\xHII)^2}   - 
\frac{2 \alpha  \xHII - \gamma  + 2 \gamma  \xHII}
{\alpha  (1-\xHII) \xHII^2 - \gamma  \xHII (1-\xHII)^2} 
\right\} d\xHII \nonumber 
\\
\tau_{\rm H} &=& 
\frac{-1}{(1-\xHII)} \Big|_{\xHIIzero}^x
- \int_{\xHIIzero}^x 
\frac{2 (\alpha+\gamma)  \xHII - \gamma  }
{\alpha  (1-\xHII) \xHII^2 - \gamma  \xHII (1-\xHII)^2} 
 d\xHII \nonumber
\\
\tau_{\rm H} &=& 
\left( \frac{1}{1-\xHIIzero} - \frac{1}{1-x} \right)
- \mathcal{I} 
\end{eqnarray}
Examining the integrand alone for a moment we have, 

\begin{eqnarray}
\frac{2 (\alpha+\gamma)  \xHII - \gamma  }
{\alpha  (1-\xHII) \xHII^2 - \gamma  \xHII (1-\xHII)^2} &=& 
% \\
\frac{2 (\alpha+\gamma)  \xHII - \gamma }
{\alpha \xHII^2 - \alpha \xHII^3 - \gamma \xHII (1 - 2\xHII + \xHII^2)} = 
\nonumber
\\
\frac{2 (\alpha+\gamma)  \xHII - \gamma }
{-(\alpha+\gamma) \xHII^3 
+ (\alpha + 2\gamma) \xHII^2 
- \gamma \xHII } &=& 
% \\
\frac{-2 \xHII + \frac{\gamma}{\alpha+\gamma} }
{\xHII^3 -\frac{\alpha + 2\gamma}{\alpha+\gamma} \xHII^2 
+\frac{\gamma}{\alpha+\gamma}\xHII} = 
\nonumber
\\
\frac{-2 \xHII + \mathcal{B} }
{\xHII \left[\xHII^2 - 
\left( 1+\mathcal{B} \right) \xHII + 
\mathcal{B} \right]} &=& 
\nonumber
\\
\frac{-2 \xHII + \mathcal{B} }
{\xHII 
\left(\xHII - 1 \right) \left( \xHII - \mathcal{B}  \right) } 
\quad &{\rm with}& \quad 
\mathcal{B} = \frac{\gamma}{\alpha+\gamma}
 = \xHII^{\rm ci,eq} < \xHII
\end{eqnarray}
Decompising this into partial fractions we have, 

\begin{eqnarray}
\left\{-2 \xHII + \mathcal{B} \right\}
\frac{1}{\xHII \left(\xHII - 1 \right) 
\left( \xHII - \mathcal{B}  \right)}
=
\left\{-2 \xHII + \mathcal{B} \right\}
\left\{
\frac{\mathcal{B}^{-1}}{\xHII}
+\frac{\left(1-\mathcal{B}\right)^{-1}}{\xHII-1}
+\frac{\left[\mathcal{B}\left(\mathcal{B}-1\right)\right]^{-1}}
{\xHII-\mathcal{B}}
\right\} = 
\nonumber
\\
\frac{-2}{\mathcal{B}} + 
\frac{1}{\xHII} + 
\frac{-2}{1-\mathcal{B}} \frac{\xHII}{\xHII-1} + 
\frac{\mathcal{B}}{1-\mathcal{B}} \frac{1}{\xHII-1} + 
\frac{-2}{\mathcal{B}(\mathcal{B}-1)} \frac{\xHII}{\xHII-\mathcal{B}} + 
\frac{1}{\mathcal{B}-1} \frac{1}{\xHII-\mathcal{B}} =
\nonumber
\\
\frac{-2}{\mathcal{B}} + 
\frac{1}{\xHII} + 
\frac{2}{1-\mathcal{B}} \frac{\xHII}{1-\xHII} + 
\frac{-\mathcal{B}}{1-\mathcal{B}} \frac{1}{1-\xHII} + 
\frac{2}{\mathcal{B}(\mathcal{B}-1)} \frac{\xHII}{\mathcal{B}-\xHII} + 
\frac{-1}{\mathcal{B}-1} \frac{1}{\mathcal{B}-\xHII} =
\end{eqnarray}
Performing the integral on each one of these terms we get, 

\begin{eqnarray}
\int_{\xHIIzero}^x \frac{-2 d\xHII}{\mathcal{B}} = 
\frac{2}{\mathcal{B}} (\xHIIzero-x)
\\
\int_{\xHIIzero}^x \frac{d\xHII}{\xHII} = \ln(x/\xHIIzero) 
\\
\frac{2}{1-\mathcal{B}} \int_{\xHIIzero}^x   
\frac{\xHII d\xHII}{1-\xHII} = 
\frac{2}{1-\mathcal{B}}
\left[ -\xHII - \ln(1-\xHII) \right] \Big|_{\xHIIzero}^{x} = 
\nonumber
\\
\frac{2}{1-\mathcal{B}}
\left[ 
(\xHIIzero - x) 
+ \ln \left( \frac{1-\xHIIzero}{1-x} \right)
\right]
\\
\frac{-\mathcal{B}}{1-\mathcal{B}} \int_{\xHIIzero}^x   
\frac{d\xHII}{1-\xHII} = 
\frac{-\mathcal{B}}{1-\mathcal{B}}
\ln \left( \frac{1-\xHIIzero}{1-x} \right)
\\
\frac{2}{\mathcal{B}(\mathcal{B}-1)} \int_{\xHIIzero}^x   
\frac{\xHII d\xHII}{\mathcal{B}-\xHII} = 
\frac{2}{\mathcal{B}(\mathcal{B}-1)} 
\left\{ 
\mathcal{B}(1-\ln \mathcal{B}^{-1}) - 
\left[ \xHII + \mathcal{B} \ln(\mathcal{B} - \xHII)
\right] \right\} 
\Big|_{\xHIIzero}^x = \nonumber
\\
\frac{2}{\mathcal{B}(\mathcal{B}-1)}
\left[ (\xHIIzero-x) + 
\mathcal{B} 
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x} 
 \right) \right]
\\
\frac{-1}{\mathcal{B}-1}
\int_{\xHIIzero}^x   
\frac{d\xHII}{\mathcal{B}-\xHII} = 
\frac{-1}{\mathcal{B}-1} 
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right)
\end{eqnarray}
The integral is then

\begin{eqnarray}
\mathcal{I} = 
\frac{2}{\mathcal{B}} (\xHIIzero-x) + 
\ln(x/\xHIIzero) +
%-------------
\frac{2}{1-\mathcal{B}}
\left[ 
(\xHIIzero - x) 
+ \ln \left( \frac{1-\xHIIzero}{1-x} \right)
\right] + 
%-------------
\nonumber
\\
\frac{-\mathcal{B}}{1-\mathcal{B}}
\ln \left( \frac{1-\xHIIzero}{1-x} \right) + 
%-------------
\frac{2}{\mathcal{B}(\mathcal{B}-1)}
\left[ (\xHIIzero-x) + 
\mathcal{B} 
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x} 
 \right) \right] + 
%--------------------
\frac{-1}{\mathcal{B}-1} 
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right)
\end{eqnarray}
Simplifying we find, 

\begin{eqnarray}
\left[
\frac{2}{\mathcal{B}} + 
\frac{2}{1-\mathcal{B}} + 
\frac{2}{\mathcal{B}(\mathcal{B}-1)}
\right] 
\left(\xHIIzero-x\right) = 
%-------------
\left[
\frac{2}{\mathcal{B}} - 
\frac{2}{\mathcal{B}-1} + 
\frac{2}{\mathcal{B}(\mathcal{B}-1)}
\right] 
\left(\xHIIzero-x\right) = 0
\\
\left[
\frac{2}{1-\mathcal{B}} + 
\frac{-\mathcal{B}}{1-\mathcal{B}}
\right]
\ln \left( \frac{1-\xHIIzero}{1-x} \right) = 
\frac{2-\mathcal{B}}{1-\mathcal{B}} 
\ln \left( \frac{1-\xHIIzero}{1-x} \right) 
\\
\left[
\frac{2 \mathcal{B}}{\mathcal{B}(\mathcal{B}-1)} - 
\frac{1}{\mathcal{B}-1}
\right]
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right) = 
\frac{1}{\mathcal{B}-1}
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right)
\end{eqnarray}
and the total inverse solution is

\begin{eqnarray}
\tau_{\rm H} = 
\left[
\frac{1}{1-\xHIIzero} - \frac{1}{1-x} 
\right]
- \ln(x/\xHIIzero)
- \frac{2-\mathcal{B}}{1-\mathcal{B}} 
\ln \left( \frac{1-\xHIIzero}{1-x} \right) 
- \frac{1}{\mathcal{B}-1}
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right)
=
\nonumber
\\
\left[
\frac{1}{1-\xHIIzero} - \frac{1}{1-x} 
\right]
+ \ln(\xHIIzero/x)
+ \frac{2-\mathcal{B}}{1-\mathcal{B}} 
\ln \left( \frac{1-x}{1-\xHIIzero} \right) 
+ \frac{1}{1-\mathcal{B}}
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right)
\end{eqnarray}
Finally, because $x$ is arbitrary, we can write,

\begin{eqnarray}
\tau_{\rm H} = 
\left[
\frac{1}{1-\xHIIzero} - \frac{1}{1-x} 
\right]
+ \frac{2-\mathcal{B}}{1-\mathcal{B}} 
\ln \left( \frac{1-x}{1-\xHIIzero} \right) 
+ \ln \left( \frac{\xHIIzero}{x} \right)
+ \frac{1}{1-\mathcal{B}}
\ln \left( 
\frac{\mathcal{B}-\xHIIzero}{\mathcal{B}-x}
\right)
\end{eqnarray}



\section{Plane Parallel Radiation on Uniform Slab with Collisions and
  Power Law Spectrum}

Here we examine the case of plane parallel radiation with a power law
spectrum incident on a semi-infinite slab of pure Hydrogen ($y = 0$)
with a constant density  $\nH$.  
In this derivation, we include the effects of collisional ionizations. 
We choose the coordinate system such that the surface of the slab is  
coincident with the $x-y$ plane and the positive $z$ direction extends
into the slab. 

Our goal is to calculate the ionization state $\xHI(z) = 1-\xHII(z)$
as a function of $z$ given a flux of incoming photons with a spectrum
of the form, 
\begin{eqnarray}
F(w) dw = F_0 w^{-\alpha} dw, \quad 
w = \left( \frac{\nu}{\nuth} \right)
\end{eqnarray}
is the number of photons per second per cm$^2$, with energies between
$w$ and $w+dw$.  
The photoionization rate is related to 
$\tau(z)$ as,
\begin{eqnarray}
\Gamma(z) = \frac{F_0 \sigma_0}{h} \int_0^q 
w^{-(4+\alpha)} e^{-\tau(z)} dw = 
\Gamma(z) = \frac{F_0 \sigma_0}{h} \int_0^q 
w^{-(4+\alpha)} e^{-\tau_0(z) w^{-3}} dw
\\
\tau(z) = \int_0^z \nH \xHI(s) \sigth w^{-3} ds = 
w^{-3} \int_0^z \nH \xHI(s) \sigth ds =
w^{-3} \tau_0(z) 
\end{eqnarray}
The fundamental theorem of calculus allows us to write, 

\begin{eqnarray}
\tau_0(z) &=& 
\int_0^z \nH \xHI(s) \sigth ds 
\\
\frac{d\tau_0(z)}{dz} &=& 
\frac{d}{dz} \int_0^z \nH \xHI(s) \sigth ds = 
\nH \xHI \sigma_0
\\
\frac{d^2\tau_0(z)}{dz^2} &=& 
\frac{d}{dz} \nH \xHI \sigma_0 = 
\nH \sigma_0 \frac{d\xHI}{dz}
\end{eqnarray}
At the surface of the slab, we have the boundary condition, 
$\tau(0) = 0$, and so,

\begin{eqnarray}
\Gamma(0) &=& \Gamma_0 = 
\frac{F_0 \sigma_0}{h} \int_1^q w^{-(4+\alpha)} dw = 
\frac{F_0 \sigma_0}{h} \left\{
\frac{w^{-(3+\alpha)}}{-(3+\alpha)}
\right\}_1^q 
\\
&=& \frac{F_0 \sigma_0}{(3+\alpha)h} \left\{
1^{-(3+\alpha)} - q^{-(3+\alpha)}
\right\} \xrightarrow{\text{q goes to $\infty$}} 
\frac{F_0 \sigma_0}{(3+\alpha)h}
\end{eqnarray}
The ionization states at $z=0$ ($\xHIzero$ and $\xHIIzero$) can be
calculated by first calculating $R$, $Q$, and $P$ and then using the solutions 
from \S 3 above with $\nH, T, \Gamma_0, y$.

\subsection{Neutral Fraction Solution}

Lets begin with the $\xHI$ solution, 

\begin{eqnarray}
\frac{d\xHI}{dt} &=& -(\Gamma + \gamma \, \nel) \xHI + 
\alpha \, \nel (1 - \xHI)
\end{eqnarray}
therefore, 

\begin{eqnarray}
(\Gamma + \gamma \nel) \xHI &=&   
\alpha \nel (1-\xHI) \nonumber 
\\
\Gamma \xHI + \gamma  \nH (1 - \xHI) \xHI  &=&   
\alpha \, \nH (1 - \xHI)^2 \nonumber 
\end{eqnarray}
therefore, 

\begin{eqnarray}
\frac{\Gamma}{\nH}  &=&   
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI}
 \nonumber 
\\
\frac{F_0 \sigma_0}{\nH h} \int_1^q 
w^{-(4+\beta)} e^{-\tau_0(z) w^{-3}} dw &=&
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI}
\nonumber
\end{eqnarray}
If we now specialize to the case of $\beta = 0$ and $q = \infty$

\begin{eqnarray}
\frac{F_0 \sigma_0}{3 \nH h \tau_0} \left( 1 - e^{-\tau_0} \right)
=
\frac{\Gamma_0}{\nH} 
\frac{\left( 1 - e^{-\tau_0} \right)}{\tau_0}   
=
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI} \nonumber 
\end{eqnarray}

\begin{eqnarray}
-e^{-\tau_0} &=& 
\left( \frac{\nH \tau_0}{\Gamma_0} \right)
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI} - 1
\nonumber 
\\
e^{-\tau_0}
&=& 1 - \left( \frac{\nH \tau_0}{\Gamma_0} \right)
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI} 
\nonumber 
\\
e^{-\tau_0}
&=& \frac{\xHI}{\xHI} - 
\frac{ (\nH/\Gamma_0) \tau_0 
\left[ 
\alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI 
\right] }
{\xHI} 
\nonumber 
\\
-\tau_0   &=& 
\ln \left\{ \frac{\xHI - 
(\nH/\Gamma_0) \tau_0 \left[ \alpha (1 - \xHI)^2 - 
                  \gamma (1 - \xHI) \xHI \right] }
{\xHI} \right\} \nonumber 
\\
-\tau_0   &=& 
\ln \left\{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\} - 
\ln \left[ \xHI \right] \nonumber 
\\
\tau_0   &=&
\ln \left[ \xHI \right] - 
\ln \left\{\xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\} 
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{d}{dz}  \ln \left[ \xHI \right] - 
\frac{d}{dz}  \ln \left\{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\} 
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\frac{d/dz \left\{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\}}
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\frac{\frac{d\xHI}{dz} 
- (\nH/\Gamma_0) \frac{d}{dz} \left\{ \tau_0 
\left[ \alpha (1 - 2\xHI + \xHI^2) 
- \gamma (\xHI - \xHI^2) \right] 
\right\}}
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\frac{\frac{d\xHI}{dz} 
- (\nH/\Gamma_0) \left( \mathcal{D}_{\alpha} - \mathcal{D}_{\gamma} \right) }
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\end{eqnarray}

%%
%% Derivative D_alpha and D_gamma
%%

\begin{eqnarray}
\mathcal{D}_{\alpha} &=& 
\frac{d}{dz}  
\left[ \alpha \tau_0 (1 - 2\xHI + \xHI^2) \right] 
\nonumber
\\
&=& 
\alpha  \left[
\frac{d\tau_0}{dz} 
- 2 \left( \frac{d\tau_0}{dz}\xHI 
+ \tau_0 \frac{d\xHI}{dz} \right) 
+ \left( \frac{d\tau_0}{dz}\xHI^2 
+ 2 \tau_0 \xHI \frac{d\xHI}{dz} \right) \right] 
\nonumber
\\
&=& 
\alpha  \left[
\left( 1 - 2\xHI + \xHI^2 \right) \frac{d\tau_0}{dz} 
+ \left( \tau_0 + 2\tau_0 \xHI \right) \frac{d\xHI}{dz} \right] 
\nonumber
\\
\mathcal{D}_{\gamma} &=& 
\frac{d}{dz}  
\left[ \gamma  \tau_0 (\xHI - \xHI^2) \right] 
\nonumber
\\
&=& 
\gamma  \left[ 
\left( \frac{d\tau_0}{dz}\xHI 
+ \tau_0 \frac{d\xHI}{dz} \right) 
- \left( \frac{d\tau_0}{dz}\xHI^2 
+ 2 \tau_0 \xHI \frac{d\xHI}{dz} \right) \right] 
\nonumber
\\
&=& 
\gamma  \left[ 
\left( \xHI - \xHI^2 \right) \frac{d\tau_0}{dz}  
+ \left( \tau_0  + 2 \tau_0 \xHI \right) \frac{d\xHI}{dz}  
\right] 
\nonumber
\\
\mathcal{D}_{\alpha} - \mathcal{D}_{\gamma} &=& 
\left[ 
\alpha  \left( 1 - 2\xHI + \xHI^2 \right)
-\gamma  (\xHI - \xHI^2) \right] \frac{d\tau_0}{dz}
+ \left( \alpha - \gamma \right)  
\left( \tau_0  + 2 \tau_0 \xHI \right) \frac{d\xHI}{dz} 
\nonumber
\end{eqnarray}


%%
%% Plug back in
%%


\begin{eqnarray}
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\nonumber
\\
&& \frac{\frac{d\xHI}{dz} 
- (\nH/\Gamma_0) \left\{ \left[ 
\alpha \left( 1 - 2\xHI + \xHI^2 \right)
-\gamma (\xHI - \xHI^2) \right] \frac{d\tau_0}{dz}
+ \left( \alpha - \gamma \right)  
\left( \tau_0  + 2 \tau_0 \xHI \right) \frac{d\xHI}{dz} \right\} }
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\nonumber
\\
&& \frac{ 
- (\nH/\Gamma_0) \left\{ \left[ 
\alpha \left( 1 - 2\xHI + \xHI^2 \right)
-\gamma (\xHI - \xHI^2) \right] \frac{d\tau_0}{dz}
+ \left[ 
\left( \alpha - \gamma \right) 
\left( \tau_0  + 2 \tau_0 \xHI \right) 
-\frac{\Gamma_0}{\nH} \right] \frac{d\xHI}{dz} \right\} }
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\end{eqnarray}
Examining the integrand alone for a moment we have, 

\begin{eqnarray}
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{\alpha (1-\xHI)^2 \xHI - \gamma (1-\xHI) \xHI^2} = 
% \\
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{\alpha (1-2\xHI+\xHI^2) \xHI - \gamma \xHI^2 + \gamma \xHI^3} = 
\nonumber
\\
\frac{2 (\alpha+\gamma) \xHI - 2 \alpha - \gamma }
{(\alpha+\gamma)\xHI^3 - (2\alpha + \gamma) \xHI^2 + \alpha \xHI} = 
%  \\
\frac{2 \xHI + \frac{-2 \alpha - \gamma}{\alpha + \gamma} }
{\xHI^3 - \frac{2\alpha + \gamma}{\alpha + \gamma} \xHI^2 + 
\frac{\alpha}{\alpha + \gamma} \xHI} = 
\nonumber
\\
\frac{2 \xHI - \frac{2 \alpha + \gamma}{\alpha + \gamma} }
{\xHI \left( \xHI^2 - \frac{2\alpha + \gamma}{\alpha + \gamma} \xHI + 
\frac{\alpha}{\alpha + \gamma} \right)} = 
% \\
\frac{2 \xHI - \left(
\frac{\alpha + \gamma}{\alpha + \gamma} + 
\frac{\alpha}{\alpha + \gamma} \right)}
{\xHI \left[ \xHI^2 - 
\left( \frac{\alpha + \gamma}{\alpha + \gamma} +
\frac{\alpha}{\alpha + \gamma} \right) 
 \xHI + 
\frac{\alpha}{\alpha + \gamma} \right]} = 
\nonumber
\\
\frac{2 \xHI - \left(
1 + \frac{\alpha}{\alpha + \gamma} \right)}
{\xHI \left[ \xHI^2 - 
\left( 1 + \frac{\alpha}{\alpha + \gamma} \right)
 \xHI + 
\frac{\alpha}{\alpha + \gamma} \right]} = 
% \\
\frac{2 \xHI - (1+\mathcal{A}) }
{\xHI \left[ \xHI^2 - (1+\mathcal{A}) \xHI + \mathcal{A} \right]} = 
\nonumber
\\
\frac{2 \xHI - (1+\mathcal{A}) }
{\xHI (\xHI - 1)(\xHI- \mathcal{A})}, 
\quad {\rm with} \quad 
\mathcal{A} = \frac{\alpha}{\alpha+\gamma}
 = \xHI^{\rm ci,eq} > \xHI
\end{eqnarray}
where $\mathcal{A}$ is the equilibrium neutral fraction in the absence
of photo-ionizations.  Now we can decompose this into partial fractions, 

\begin{eqnarray}
\frac{2 \xHI - (1+\mathcal{A}) }
{\xHI (\xHI - 1)(\xHI- \mathcal{A})} &=& 
\nonumber
\\
\left\{ 2 \xHI - (1+\mathcal{A}) \right\}
\left\{ 
\frac{\mathcal{A}^{-1}}{\xHI} + 
\frac{(1-\mathcal{A})^{-1}}{\xHI-1} + 
\frac{[\mathcal{A}(\mathcal{A}-1)]^{-1}}{\xHI-\mathcal{A}} 
\right\} &=& 
\nonumber
\\
\frac{2\mathcal{A}^{-1}}{1} -
\frac{\mathcal{A}^{-1} + 1}{\xHI} + 
\frac{2\xHI(1-\mathcal{A})^{-1}}{\xHI-1} - 
\frac{(1+\mathcal{A})/(1-\mathcal{A})}{\xHI-1} + 
\frac{2\xHI[\mathcal{A}(\mathcal{A}-1)]^{-1}}{\xHI-\mathcal{A}} - 
\frac{(1+\mathcal{A})[\mathcal{A}(\mathcal{A}-1)]^{-1}}{\xHI-\mathcal{A}} 
&=& \nonumber
\\
\frac{2}{\mathcal{A}} -
\frac{\mathcal{A}^{-1} + 1}{\xHI} - 
\frac{2}{1-\mathcal{A}} \frac{\xHI}{1-\xHI} + 
\frac{1+\mathcal{A}}{1-\mathcal{A}} \frac{1}{1-\xHI} - 
\frac{2}{\mathcal{A}(\mathcal{A}-1)} \frac{\xHI}{\mathcal{A}-\xHI} + 
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\frac{1}{\mathcal{A}-\xHI} 
&&
\end{eqnarray}
In the last step I've swapped negative signs on some of the terms so
that the individual integrations yield non-imaginary values.  Now, we can
integrate these term by term,

\begin{eqnarray}
\frac{2}{\mathcal{A}} \int_{\xHIzero}^x d\xHI  = 
\frac{2}{\mathcal{A}} (x - \xHIzero)
\\
-(\mathcal{A}^{-1} + 1)\int_{\xHIzero}^x \frac{d\xHI}{\xHI} = 
-(\mathcal{A}^{-1} + 1) \ln({x/\xHIzero})
\\
\frac{-2}{1-\mathcal{A}} 
\int_{\xHIzero}^x 
\frac{\xHI d\xHI}{1-\xHI}  = 
\frac{-2}{1-\mathcal{A}} 
[ -\xHI - \ln(1-\xHI) ] \Big|_{\xHIzero}^x = \nonumber
\\
\frac{-2}{1-\mathcal{A}} 
\left[ (\xHIzero-x) + \ln \left( \frac{1-\xHIzero}{1-x} \right) \right]
\\
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\int_{\xHIzero}^x \frac{d\xHI}{1-\xHI} = 
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\ln \left( \frac{1-\xHIzero}{1-x} \right) 
\\
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\int_{\xHIzero}^x \frac{\xHI d\xHI}{\mathcal{A}-\xHI} = 
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\left\{ 
\mathcal{A}(1-\ln \mathcal{A}^{-1}) - 
\left[ \xHI + \mathcal{A} \ln(\mathcal{A} - \xHI)
\right] \right\} 
\Big|_{\xHIzero}^x = \nonumber
\\
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\left[ (\xHIzero-x) + 
\mathcal{A} 
\ln \left( 
\frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} 
 \right) \right]
\\
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\int_{\xHIzero}^x \frac{d\xHI}{\mathcal{A}-\xHI} = 
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\ln \left( 
\frac{\mathcal{A} - \xHIzero}{\mathcal{A}-x}
\right)
\end{eqnarray}
The integral $\mathcal{I}$ is then,

\begin{eqnarray}
\mathcal{I} = 
\frac{2}{\mathcal{A}} (x - \xHIzero) - 
(\mathcal{A}^{-1} + 1) \ln({x/\xHIzero}) + 
\nonumber
\\
\frac{-2}{1-\mathcal{A}} 
\left[ (\xHIzero-x) + 
\ln \left( \frac{1-\xHIzero}{1-x} \right) \right] + \nonumber
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\ln \left( \frac{1-\xHIzero}{1-x} \right) +
\\
\frac{-2}{\mathcal{A}(\mathcal{A}-1)}
\left[ (\xHIzero-x) + 
\mathcal{A} 
\ln \left( 
\frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} 
 \right) \right] + 
\frac{1+\mathcal{A}}{\mathcal{A}(\mathcal{A}-1)}
\ln \left( 
\frac{\mathcal{A} - \xHIzero}{\mathcal{A}-x}
\right) 
\end{eqnarray}
Rearranging negative signs such that terms in $()$ are positive,

\begin{eqnarray}
\mathcal{I} = \frac{2}{\mathcal{A}} (x - \xHIzero) - 
(1 + \mathcal{A}^{-1}) \ln({x/\xHIzero}) + 
\nonumber
\\
\frac{2}{1-\mathcal{A}} 
\left[ (x - \xHIzero) - 
\ln \left( \frac{1-\xHIzero}{1-x} \right) \right] + \nonumber
\frac{1+\mathcal{A}}{1-\mathcal{A}} 
\ln \left( \frac{1-\xHIzero}{1-x} \right) +
\\
\frac{-2}{\mathcal{A}(1-\mathcal{A})}
\left[ (x-\xHIzero) - 
\mathcal{A} 
\ln \left( 
\frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} 
 \right) \right] - 
\frac{1+\mathcal{A}}{\mathcal{A}(1-\mathcal{A})}
\ln \left( 
\frac{\mathcal{A} - \xHIzero}{\mathcal{A}-x}
\right) = 
\end{eqnarray}
Simplifying we see,

\begin{eqnarray}
\left[ 
\frac{2}{\mathcal{A}} + 
\frac{2}{1-\mathcal{A}} -
\frac{2}{\mathcal{A}(1-\mathcal{A})} 
\right] (x-\xHIzero) = 0
\\
\left[ 
-\frac{2}{1-\mathcal{A}} + 
\frac{1+\mathcal{A}}{1-\mathcal{A}}
\right] 
\ln \left( \frac{1-\xHIzero}{1-x} \right) = 
%-----------------
- 
\ln \left( \frac{1-\xHIzero}{1-x} \right)
\\
\left[ 
\frac{2\mathcal{A}}{\mathcal{A}(1-\mathcal{A})} -
\frac{1+\mathcal{A}}{\mathcal{A}(1-\mathcal{A})} 
\right] 
\ln \left( \frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} \right) = 
%-----------------------------
\left[ 
-\frac{1}{\mathcal{A}} 
\right] 
\ln \left( \frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} \right) 
\end{eqnarray}
and the total inverse solution is, 

\begin{eqnarray}
\tau_{\rm H} &=& 
\left( \frac{1}{\xHIzero} - \frac{1}{x} \right) +
(1 + \mathcal{A}^{-1}) 
\ln \left( \frac{x}{\xHIzero} \right) +
\ln \left( \frac{1-\xHIzero}{1-x} \right) + 
\frac{1}{\mathcal{A}} 
\ln \left( \frac{\mathcal{A}-\xHIzero}{\mathcal{A}-x} \right) 
\\
&=&
\left( \frac{1}{\xHIzero} - \frac{1}{x} \right) +
\ln \left[ \frac{x(1-\xHIIzero)}{\xHIzero(1-x)} \right] + 
\frac{1}{\mathcal{A}}
\ln \left[ \frac{x(\mathcal{A}-\xHIIzero)}{\xHIzero(\mathcal{A}-x)} \right] 
\end{eqnarray}
Finally, because $x$ is arbitrary, we can write,

\begin{eqnarray}
\tau_{\rm H} = 
\left( \frac{1}{\xHIzero} - \frac{1}{\xHI} \right) +
\ln \left[ \frac{\xHI(1-\xHIIzero)}{\xHIzero(1-\xHI)} \right] + 
\frac{1}{\mathcal{A}}
\ln \left[ 
\frac{\xHI(\mathcal{A}-\xHIIzero)}
{\xHIzero(\mathcal{A}-\xHI)} 
\right] 
\end{eqnarray}



\section{Plane Parallel Radiation on Uniform Slab with Collisions and
  Arbitrary Spectrum - Tabular Solution}

Here we examine the case of plane parallel radiation with an arbitrary
spectrum incident on a semi-infinite slab of pure Hydrogen ($y = 0$)
with a constant density  $\nH$.  In this derivation, we include the
effects of collisional ionizations.  We choose the coordinate system
such that the surface of the slab is coincident with the $x-y$ plane
and the positive $z$ direction extends into the slab.  


\subsection{Neutral Fraction Solution}

Lets begin with the $\xHI$ equilibrium solution, 

\begin{eqnarray}
\frac{d\xHI}{dt} = 0 = -(\Gamma + \gamma \, \nel) \xHI + 
\alpha \, \nel (1 - \xHI) 
\end{eqnarray}
therefore, 

\begin{eqnarray}
(\Gamma + \gamma \nel) \xHI &=&   
\alpha \nel (1-\xHI) \nonumber 
\\
\Gamma \xHI + \gamma  \nH (1 - \xHI) \xHI  &=&   
\alpha \, \nH (1 - \xHI)^2 \nonumber 
\\
\frac{\Gamma}{\nH}  &=&   
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI}
 \nonumber 
\end{eqnarray}
The photoionization rate can be written, 

\begin{eqnarray}
\Gamma  &=&   
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI}
 \nonumber 
\\
\frac{F_0 \sigma_0}{\nH h} \int_1^q 
w^{-(4+\beta)} e^{-\tau_0(z) w^{-3}} dw &=&
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI}
\nonumber
\end{eqnarray}
If we now specialize to the case of $\beta = 0$ and $q = \infty$

\begin{eqnarray}
\frac{F_0 \sigma_0}{3 \nH h \tau_0} \left( 1 - e^{-\tau_0} \right)
=
\frac{\Gamma_0}{\nH} 
\frac{\left( 1 - e^{-\tau_0} \right)}{\tau_0}   
=
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI} \nonumber 
\end{eqnarray}

\begin{eqnarray}
-e^{-\tau_0} &=& 
\left( \frac{\nH \tau_0}{\Gamma_0} \right)
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI} - 1
\nonumber 
\\
e^{-\tau_0}
&=& 1 - \left( \frac{\nH \tau_0}{\Gamma_0} \right)
\frac{\alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI }
{\xHI} 
\nonumber 
\\
e^{-\tau_0}
&=& \frac{\xHI}{\xHI} - 
\frac{ (\nH/\Gamma_0) \tau_0 
\left[ 
\alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI 
\right] }
{\xHI} 
\nonumber 
\\
-\tau_0   &=& 
\ln \left\{ \frac{\xHI - 
(\nH/\Gamma_0) \tau_0 \left[ \alpha (1 - \xHI)^2 - 
                  \gamma (1 - \xHI) \xHI \right] }
{\xHI} \right\} \nonumber 
\\
-\tau_0   &=& 
\ln \left\{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\} - 
\ln \left[ \xHI \right] \nonumber 
\\
\tau_0   &=&
\ln \left[ \xHI \right] - 
\ln \left\{\xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\} 
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{d}{dz}  \ln \left[ \xHI \right] - 
\frac{d}{dz}  \ln \left\{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\} 
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\frac{d/dz \left\{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - \gamma (1 - \xHI) \xHI \right] 
\right\}}
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\frac{\frac{d\xHI}{dz} 
- (\nH/\Gamma_0) \frac{d}{dz} \left\{ \tau_0 
\left[ \alpha (1 - 2\xHI + \xHI^2) 
- \gamma (\xHI - \xHI^2) \right] 
\right\}}
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\\
\frac{d \tau_0}{dz}   &=& 
\frac{1}{\xHI} \frac{d\xHI}{dz} - 
\frac{\frac{d\xHI}{dz} 
- (\nH/\Gamma_0) \left( \mathcal{D}_{\alpha} - \mathcal{D}_{\gamma} \right) }
{ \xHI - (\nH/\Gamma_0) \tau_0 
\left[ \alpha (1 - \xHI)^2 - 
\gamma (1 - \xHI) \xHI \right] }  
\nonumber 
\end{eqnarray}





\section{Stromgren Sphere}

This classic problem is that of a monochromatic point source of
radiation emitting isotropically into a uniform medium of pure Hydrogen.  
If we choose the temperature of this gas to be low enough we can set
both $\gamma$ and $y$ to zero.  We choose the coordinate system such that
the point source is at the origin.  Our goal is to calculate
the ionization fraction $\xHII(r)$ as a function of $r$ given a flux
of outgoing photons per second, $F_0$.  In order to accomplish this we
must calculate the optical depth at the hydrogen ionizing threshold, 
$\tau(r)$, as well. 

The photoionization rate is related to 
$\tau(r) = \int_0^r \nH (1-\xHII) \sigth ds$ as, 
\begin{eqnarray}
\Gamma(r) = F(r) \sigth = 
\frac{F_0}{4 \pi r^{2}} e^{-\tau(r)}\sigth 
\end{eqnarray} 
At the position of the source we have infinite flux, but we can choose
some small $r = r_0$ to act as our inner most radius.  At this radius
we have a boundary condition, $\tau(r=r_0) = 0$, and so 
$\Gamma(r=r_0) = \Gamma_0 = F_0 \sigth / (4 \pi r_0^2)$.
With $\gamma$ and $y$ equal to zero, the equilibrium ionization
fraction obeys the equation,
\begin{eqnarray}
\frac{\xHII^2}{1-\xHII} = \frac{\Gamma}{\alpha \nH} = 
\frac{F_0 \sigth {\rm Exp}[{-\int_0^r \nH (1-\xHII) \sigth ds}]}
{\alpha \nH 4 \pi r^2}
\end{eqnarray}
First let us solve for the ionization fraction at $r=r_0$ or $x_0$,
\begin{eqnarray}
\frac{x_0^2}{1-x_0} = 
\frac{F_0 \sigth}{\alpha \nH 4 \pi r_0^2} = \mathcal{C} > 0, \quad
x_0^2 + \mathcal{C} x_0 - \mathcal{C} = 0
\end{eqnarray}
Taking the positive root of the quadratic formula keeps $x_0$ positive
and so we have,
\begin{eqnarray}
x_0 = \frac{-\mathcal{C} + \sqrt{\mathcal{C}^2 + 4\mathcal{C}}}{2}
\end{eqnarray}


Now lets move on to solving for the $z$ dependent ionization
fraction.  We have,
\begin{eqnarray}
\frac{\xHII^2}{1-\xHII} \frac{r^2}{\mathcal{C}} \frac{1}{r_0^2}&=& 
{\rm Exp} \left[ \int_{r_0}^r \nH (1-\xHII) \sigth ds \right] \\
\ln \left[ \frac{\xHII^2}{1-\xHII} \frac{r^2}{\mathcal{C}} 
\frac{1}{r_0^2} \right] &=& 
\int_{r_0}^r \nH (1-\xHII) \sigth ds \\
\frac{d}{dr} \left[ \ln(\xHII^2) - \ln(1-\xHII) - 
\ln(r_0^2\mathcal{C}) + \ln(r^2) \right]&=& 
\frac{d}{dr} \int_{r_0}^r \nH (1-\xHII) \sigth ds \\
\frac{d}{dr} \left[ 2\ln(\xHII) - \ln(1-\xHII) + 2\ln(r) \right]&=& 
\nH (1-\xHII) \sigth \\
\frac{2}{\xHII} \frac{d\xHII}{dr} + \frac{1}{1-\xHII}
\frac{d\xHII}{dr} + \frac{2}{r} &=&
\nH (1-\xHII) \sigth \\
\left[ \frac{2}{\xHII (1-\xHII)} + \frac{1}{(1-\xHII)^2} \right]
\frac{d\xHII}{dr} + \frac{2}{r(1-\xHII)}&=& \nH \sigth \\
f^{-1}(\xHII) \frac{d\xHII}{dr} + \frac{2}{r(1-\xHII)}
&=& \nH \sigth \\
\frac{d\xHII}{dr} + \frac{2f(\xHII)}{r(1-\xHII)} - \nH \sigth f(\xHII)
&=& 0 
\end{eqnarray}

This relation gives the distance from the surface of the slab in terms
of the ionized fraction $x = \xHII$.  This is the inverse of the
function we were looking for, but it can be used to map out the solution.
Starting from $x=x_0$ we can decrease x in as many increments as we
need to cover a suitable range in $z$. 







\bibliographystyle{alpha} % or "unsrt", "alpha", "abbrv", etc.
\bibliography{biblio}	  % use data in file "astrobibl.bib"


\end{document}
